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The  Fock-space  coupled  cluster  (FSCC)  approach  to  the  electron  correlation  problem 
offers  several  advantages  over  the  more  familiar  single  reference  coupled  cluster  (CC) 
methods.  It  is  capable  of  treating  states  with  a  multireference  character,  and  it  is  capable 
of  providing  state-to-state  energy  differences  for  several  states  in  a  single  calculation.  As 
explained  in  this  work,  basic  differences  between  the  two  theories  result  in  higher-excited 
cluster  operators  being  more  important  in  FSCC  methods  than  CC  methods. 

In  order  to  improve  the  accuracy  of  FSCC  calculations,  this  work  involves  the 
addition  of  triple  excitation  effects  to  the  widely  applied  FSCC  model  restricted  to 
single  and  double  excitations  (FSCCSD)  for  ionization  potentials.  Equations  for  the 
full  model,  including  single,  double,  and  triple  excitations  (FSCCSDT),  are  derived  and 
various  approximations  are  proposed  as  a  compromise  between  the  improvement  in  the 
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wavefunction  and  energy  due  to  the  inclusion  of  triples  and  the  computational  costs, 
which  are  significantly  larger  than  for  the  FSCCSD  model. 

Two  approximations  to  the  full  FSCCSDT  model  have  been  implemented  and  their 
performance  is  examined  for  a  variety  of  molecules.  Results  are  found  to  be  generally 
good,  providing  results  comparable  to  single  reference  coupled  cluster  methods  with 
approximate  inclusion  of  triple  excitations. 

Finally,  both  CC  and  FSCC  methods  are  used  to  study  the  structures  of  several 
small  carbon  clusters,  a  class  of  molecules  which  has  been  shown  to  require  high-quality 
theoretical  methods  and  which  can  make  use  of  both  FSCC's  multireference  and  direct 
energy  difference  capabilities.  Detailed  calculations  show  that  the  ground  state  of  Cs+  is  a 
2E+,  with  Cqov  symmetry,  but  with  low  frequency  bending  vibrational  modes.  Additional 
calculations  for  other  odd  carbon  cluster  cations  suggest  that  the  "linear"  structures  of 
C7+  and  C9+  observed  experimentally  may  in  fact  be  distorted  from  true  linearity  due 
to  the  Renner-Teller  effect. 
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CHAPTER  1 
INTRODUCTION 


In  recent  years,  advances  in  both  experimental  and  theoretical  chemistry  have  fueled 
an  exciting  synergy  between  the  two  fields  [1].  Experimental  innovations  have  allowed 
study  of  molecules  in  conditions  similar  to  those  commonly  employed  in  theoretical 
calculations  (low  temperature,  low  density),  producing  data  of  increasing  accuracy  and 
specificity.  Likewise,  through  increases  in  computational  power  and  methodological 
sophistication,  theoretical  chemists  have  been  able  to  tackle  increasingly  complex  systems 
with  high  accuracy.  Where  experimental  and  theoretical  data  are  direcdy  comparable, 
critical  evaluations  of  both  approaches  are  possible.  In  other  cases,  the  data  are  often 
complementary,  allowing  one  or  the  other  researcher  to  bring  the  problem  into  finer  focus. 

This  has  been  especially  apparent  in  research  on  atomic  clusters,  a  field  which  has 
expanded  rapidly  in  the  past  decade.  Clusters  are  interesting  for  a  variety  of  rea- 
sons. Metal  clusters  are  used  as  models  for  the  study  of  catalytic  processes.  Clus- 
ters of  lighter  atoms,  particularly  carbon,  have  been  identified  in  circumstellar  space 
and  are  expected  also  to  be  found  in  interstellar  molecular  clouds,  where  they  may 
be  important  in  the  creation  of  new  stars.  They  have  also  been  found  in  combustion 
products  on  Earth.  Clusters  may  have  useful  and  unique  physical  properties,  mak- 
ing them  interesting  from  a  materials  viewpoint,  and  they  are  studied  as  models  of 
how  properties  of  bulk  materials  are  approached  by  increasingly  large  collections  of 
atoms. 
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2 
Although  the  literature  on  clusters  is  expanding  rapidly,  the  sought  after  understanding 
has  not  yet  emerged.  One  reason  for  this  is  that  clusters,  even  relatively  small  ones,  are 
rather  unlike  "traditional"  molecules  which  scientists  have  studied  for  centuries,  and 
they  often  defy  the  "rules"  and  intuitions  which  researchers  try  to  bring  from  traditional 
chemistry  to  cluster  chemistry.  Consequently,  workers  in  the  area  are  having  to  develop 
new  ideas  and  new  "rules"  about  the  structure  and  properties  of  clusters. 

An  important  component  of  creating  this  new  understanding  is  having  reliable  data 
on  the  structure  and  properties  of  many  clusters.  Since  there  are  few  experiments  which 
give  direct,  unambiguous  structural  data,  theoretical  approaches  are  particularly  useful  in 
this  context.  When  the  capabilities  of  experiment  and  theory  do  overlap,  it  is  possible 
to  compare  results  critically  and  perhaps  improve  both.  To  be  truly  useful,  however, 
the  theoretical  methods  employed  must  be  accurate  and  reliable.  The  subject  of  this 
dissertation,  the  inclusion  of  triple  excitation  effects  in  the  Fock-space  coupled  cluster 
(FSCC)  method,  was  undertaken  in  response  to  exactly  this  need. 

Coupled  cluster  (CC)  and  related  many-body  perturbation  theory  (MBPT)  methods 
(see  Chapter  3)  have  proven  their  ability  to  deliver  extremely  high-quality  properties  and 
energetics  for  a  wide  range  of  complex  systems,  and  because  of  their  size  extensivity 
(proper  scaling  of  the  energy  with  the  size  of  the  system  under  study)  [2,3],  are  preferred 
over  configuration  interaction  (CI)  methods  for  theoretical  studies  of  clusters.  Not 
every  situation,  however,  is  ideally  suited  to  the  use  of  traditional  single  reference  CC 
methods.  For  example,  close-lying  electronic  states  may  cause  problems  for  even  the  most 
sophisticated  CC  methods,  thus  requiring  a  multireference  description  of  the  wavefunction 


3 
to  obtain  meaningful  results.  Ionization  potentials  (IPs)  and  electron  affinities  (EAs)  are 
often  obtained  using  experimental  studies  of  clusters  in  photoionization  and  photoelectron 
spectroscopies.  In  order  to  provide  comparable  theoretical  results  for  interpretation  or 
confirmation  of  the  spectra,  separate  calculations  for  the  ground  state  and  each  excited 
state  of  interest  are  required  when  using  familiar  CC  methods.  It  would  be  useful  to 
have  a  method  capable  of  giving  IPs  and  EAs  for  several  states  of  interest  in  a  single 
calculation.  Fock-space  coupled  cluster  methods  are  one  way  to  approach  both  of  these 
problems. 

Often  referred  to  as  Fock-space  multireference  coupled  cluster  theory,  FSCC  (see 
Chapter  4)  can  provide  a  multireference  description  of  chosen  electronic  states  of  a 
molecule.  Considered  in  another  context,  it  can  calculate  the  energies  of  several  electron- 
detached  or  -attached  states  at  the  same  time,  thus  simplifying  the  calculation  of  IPs  and 
EAs.  Because  of  formal  differences  between  single  reference  CC  and  FSCC  theories, 
however,  a  given  truncation  of  the  two  methods  does  not  necessarily  give  results  of 
equivalent  quality.  For  example,  limiting  both  methods  to  single  and  double  excitations, 
we  find  that  an  FSCCSD  description  of  a  particular  state  may  leave  out  contributions 
which  could  be  present  in  the  CCSD  description  of  the  same  state.  Consequently,  in 
order  to  improve  the  quality  of  the  FSCC  description,  this  work  examines  the  inclusion 
of  triple  excitation  effects  in  the  FSCC  method. 

After  setting  the  stage  with  a  review  of  theory  and  terminology  leading  up  to  the 
correlation  problem,  the  traditional  coupled  cluster  equations  are  presented,  emphasizing 
their  use  as  an  "input"  to  the  FSCC  method.  Then,  in  some  detail,  equations  for  the  FSCC 
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method  are  derived,  focusing  on  the  FSCCSDT  model  for  IPs.  The  importance  of  triple 
excitations  is  examined,  in  comparison  with  single  reference  CCSD  and  CCSDT  methods, 
and  several  approximations  are  proposed  as  a  compromise  between  the  computational  cost 
of  evaluating  the  complete  triples  contribution  and  quality  required  of  the  wavefunction 
and  energies.  Although  this  development  was  spurred  by  an  interest  in  carbon  clusters, 
FSCC  methods,  with  or  without  triple  excitations,  are  by  no  means  limited  to  this 
application.  In  order  to  assess  the  performance  of  the  FSCC  methods,  calculations  are 
presented  for  a  variety  of  molecules.  Finally,  Fock-space  and  single  reference  coupled 
cluster  methods  are  used  to  study  several  carbon  cluster  cations. 


CHAPTER  2 
THEORETICAL  BACKGROUND 


2.1  Overview 

Much  of  the  theoretical  chemistry  literature  assumes  the  reader  is  a  practitioner  in  the 
field.  This  chapter  is  intended  to  provide  a  common  ground  for  the  concepts  and  notation 
used  throughout  this  dissertation.  It  is  by  no  means  exhaustive,  but  should  help  non- 
specialists  to  understand  the  material  that  follows.  References  are  given  to  the  literature 
for  more  detailed  treatments  of  various  topics. 

2.2  The  Schrodinger  Equation 

We  seek  to  solve  the  time-independent  Schrodinger  equation, 

ff|#)  =  £|tf)  (2.1) 

for  a  molecular  system  (one  or  more  molecules  of  interest  isolated  from  interaction  with 
anything  else).  This  equation  tells  us  that  a  system  in  an  eigenstate  |^)  of  the  Hamiltonian 
operator  H  has  an  energy  E.  The  usual  nonrelativistic  Hamiltonian  can  be  divided  into 
five  different  terms, 

H  =  Te  +  Tn  +  Ven  +  Vee  +  Vnn  (2.2) 

representing,  respectively,  the  electronic  and  nuclear  kinetic  energies  and  the  potentials 
due  to  interaction  of  electrons  and  nuclei,  electrons  with  other  electrons,  and  nuclei  with 
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other  nuclei.  The  individual  terms,  in  Hartree  atomic  units,  are  defined  as 
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I 


2MA 


f..M)  =  -EEt' 


i       A 


Ra\  <2-3> 


"-M-iEtnr 


"-<*>  4  E  ^% 


2  ^— '  \R a  —  Rr\ 

In  these  expressions,  i  and  j  refer  to  electrons,  A  and  5  refer  to  nuclei.  Ra  and  r,  refer  to 
the  spatial  coordinates  of  nucleus  A  and  electron  i,  respectively,  and  ZA  and  Ma  are  the 
charge  and  mass  of  nucleus  A.  The  unsubscripted  symbols  r  and  R  refer  to  the  complete 
set  of  position  vectors  of  the  electrons  and  nuclei,  respectively. 

2.3  Potential  Energy  Surfaces 

The  terms  of  the  full  Hamiltonian,  Eq.    2.2,  can  be  grouped  according  to  their 
dependence  on  the  positions  of  the  nuclei  into  a  purely  electronic  component, 

He  =  Te  +  Vee,  (2.4) 

a  purely  nuclear  term, 

Hn  =  Tn,  (2.5) 
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and  a  potential  which  couples  them, 

Vint(r,  R)  =  Ven(r,  R)  +  Vnn(R).  (2.6) 

If  we  choose  one  set  of  positions  R  as  the  origin,  denoted  by  0,  we  can  expand  the 
electron-nuclear  interaction  potential  near  this  point  in  a  Taylor  series, 

K„«(^,  =  K„,(,0)  +  X:(^)^  +  iE(^)o^  +  ....    (2.7, 

The  first  term  in  the  expansion  can  be  thought  of  as  the  potential  felt  by  the  electrons 
due  to  a  fixed  set  of  nuclei,  leading  to  the  electronic  Schrodinger  equation, 

(Re(r)  +  ViB«(r,0))|*jfc(r))  -  E'k\$k(r)),  (2.8) 

for  a  particular  state  \®k(r))  Wlth  energy  E'k. 

In  order  to  see  how  varying  the  nuclear  positions  effects  these  solutions,  we  must 
return  to  the  full  Schrodinger  equation  2.1.  By  expanding  the  full  wavefunction  in  terms 
of  the  solutions  to  the  electronic  wavefunction  and  a  nuclear  component, 

|*(r,i2))=£|*fc(r))Xfc(/2),  (2.9) 

k 

and  recognizing  that  the  Hamiltonian  of  Eq.    2.8  is  a  subset  of  the  full  Hamiltonian, 
we  obtain 

{He  +  Vintfi  +  Hn  +  W)  Y,  \$k)Xk  =  EY,  l^>X*'  (110) 

k  k 

where  V^^u  denotes  the  leading  term  in  the  Vini  expansion,  and  W denotes  the  remainder, 
W(r,  R)  —  V,nt{r,  R)  —  Vint(r,  0),  known  as  the  vibronic  interaction.  He  +  Vtntto  is  often 
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referred  to  as  the  electrostatic  Hamiltonian,  and  is  dependent  on  the  nuclear  coordinates 
only  parametrically.    Projecting  on  the  left  with  an  electronic  eigenfunction,  ($m(r)|, 
we  obtain  a  coupled  set  of  equations,  (note  that  Hn  is  independent  of  the  electronic 
coordinates  r), 

(E'm  +  Hn)Xm  +  X)(*m|W|*ft>Xt  =  EXm,  (2.11) 

k 

or  using  the  more  compact  matrix  element  notation  for  W, 

{E  -  E'm  -  Hn)Xm  =  J2  WmkXk.  (2.12) 

k^m 

If  there  is  no  vibronic  interaction  between  electronic  states  (Wmk  =  0),  then  the 
equations  decouple,  giving 

(Hn(R)  +  E'm)Xm(R)  =  EXm(R)  V  m,  (2.13) 

which  is  the  Schrodinger  equation  for  the  nuclei  moving  in  the  mean  field  of  electronic 
state  |$m(r)),  and  the  solution  to  the  full  Schrodinger  equation  is  simply  \ty(r,R))  = 
|$m(r))Xm  (/?.).  This  is  the  origin  of  the  idea  of  a  potential  energy  surface  (PES)  as  the  set 
of  nuclear  configurations  which  satisfy  Eq.  2.13  for  a  particular  electronic  state.  Taking 
Wmk  =  0  is  known  as  the  simple  adiabatic  or  Born-Oppenheimer  approximation  [4,5]. 

In  practice,  although  Wmk  is  not  in  general  zero,  the  Born-Oppenheimer  approxima- 
tion is  usually  quite  good  because  of  the  fact  that  the  nuclei  of  a  molecule  are  about  1836 
times  more  massive  than  the  electrons,  so  we  can  usually  think  of  the  nuclei  moving 
slowly  in  the  average  field  of  the  electrons,  which  are  able  to  adapt  almost  instanta- 
neously to  the  nuclear  motions.  This  approximation  is  used  in  the  majority  of  electronic 
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structure  calculations  to  date,  which  simply  evaluate  the  electronic  wavefunction  and 
energy,  Eq.  2.8,  at  various  points  on  the  PES  (nuclear  configurations)  and  ignore  the 
nuclear  dynamics. 

2.4  Calculations  on  Potential  Energy  Surfaces 

The  Born-Oppenheimer  potential  energy  surface  for  a  particular  state  of  a  molecule 
with  N  nuclei  is  (neglecting  invariance  to  rigid  translations  and  rotations  of  the  entire 
nuclear  framework)  a  .^-dimensional  space  which  can  be  characterized  according  to  its 
local  shape.  A  stationary  point  is  one  at  which  the  gradient  of  the  energy  with  respect  to 
nuclear  displacements  be  zero.  Stationary  points  may  be  minima,  transition  states,  higher- 
order  saddle  points,  or  maxima.  If  the  curvature  (second  derivative)  at  a  stationary  point 
is  positive,  the  point  is  a  local  minimum,  which  is  a  locally  stable  geometric  conformation 
for  the  state  in  question.  If  the  curvature  in  one  of  the  3N  directions  is  negative,  with 
all  others  positive,  the  point  is  a  transition  state,  which  is  a  local  maximum  on  a  path 
connecting  two  local  minima;  transition  states  are  very  important  in  the  theory  of  chemical 
reactions.  If  the  stationary  point  has  more  than  one  direction  of  negative  curvature,  it  is 
a  saddle  point  (or  local  maximum),  which  is  not  chemically  meaningful. 

Stationary  points  on  potential  energy  surfaces  can  be  located  using  numerical  opti- 
mization techniques,  such  as  Newton-Raphson,  quasi-Newton-Raphson,  etc.,  but  because 
individual  calculations  of  the  energy,  gradient,  and  hessian  (matrix  of  second  derivatives 
of  the  energy),  as  required  by  the  optimization  method,  are  extremely  expensive  relative 
to  evaluation  of  the  optimization  step,  it  is  usually  necessary  to  use  chemical  intuition 
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to  choose  regions  of  the  PES  on  which  to  focus  and  sometimes  to  restrict  the  search  to 
a  subset  of  the  full  ^-dimensional  space. 

The  hessian,  which  allows  a  stationary  point  to  be  characterized  as  a  minimum, 
transition  state,  or  saddle  point,  is  also  the  primary  input  for  evaluation  of  the  vibrational 
frequencies  at  that  structure  (in  the  harmonic  approximation).  Since  different  molecules 
and  different  structures  will  generally  have  different  vibrational  spectra,  it  serves  as  a 
useful  tool  for  the  identification  and  characterization  of  unknown  species. 

If  the  electronic  state  of  the  molecule  of  interest  is  not  known  from  other  information, 
it  may  be  necessary  to  consider  potential  energy  surfaces  for  several  states  in  order  to 
locate  the  lowest  energy  global  minimum  of  all  the  states,  which  is  known  as  the  ground 
state  of  the  system. 

2.5  Vibronic  Interactions 


Although  the  Born-Oppenheimer  approximation  and  the  idea  of  potential  energy 
surfaces  have  served  quantum  chemistry  well,  there  are  certainly  cases  in  which  the  Wml 
matrix  elements  connecting  different  electronic  states  are  significantly  different  from  zero, 
and  the  Bom-Oppenheimer  picture  is  no  longer  adequate. 

Consider  the  vibronic  interaction  term,  W,  truncated  after  the  quadratic  term,  which 
is  adequate  for  purposes  of  explanation.  Making  a  linear  transformation  from  cartesian 
coordinates   R  to  symmetrized  normal  coordinates  Q  (obtained  by  diagonalizing  the 
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hessian  matrix),  we  have 


Each  normal  coordinate,  QQ,  belongs  to  an  irreducible  representation  (irrep)  TQ  of  the 
molecular  point.  Likewise,  the  electronic  states  |$m)  belong  to  an  irrep,  Tm.  Thus  from 
group  theory,  we  know  that  the  linear  vibronic  constant, 

'dVint 


$m>,  (2-15) 


is  nonzero  if  and  only  if  Ta  e  [I\-  x  Tm].  This  is  a  general  result  for  the  interaction 
of  any  two  states,  but  of  particular  importance  is  the  case  in  which  the  two  electronic 
states  belong  to  the  same  irreducible  representation.  If  the  irrep  is  not  degenerate,  then 
[Tk  x  Tk]  =  I\-  x  I\.  =  r^j,  so  that  the  coupling  is  via  totally  symmetric  motions,  which 
by  definition  cannot  change  the  symmetry  of  the  molecule.  If,  on  the  other  hand,  the 
two  electronic  states  (or  components  of  states)  belong  to  a  degenerate  irrep,  the  direct 
product  results  in  non-totally  symmetric  representations,  and  thus  the  states  are  coupled 
by  motions  which  break  (lower)  the  symmetry  of  the  molecular  framework.  This  is  the 
basis  of  the  Jahn-Teller  (JT)  theorem  [6],  which  says  that  for  a  nonlinear  polyatomic 
molecule  in  an  n-fold  degenerate  electronic  state,  the  high  symmetry  point  at  which  the 
n  potential  surfaces  intersect  cannot  be  a  minimum  on  all  n  surfaces.  In  other  words, 
there  will  always  be  at  least  one  direction  for  which  the  energy  is  lowered  when  the 
degeneracy  is  lifted  (by  moving  away  from  the  high  symmetry  point).  The  tendency  of 
a  molecule  in  a  degenerate  state  to  lower  its  symmetry  and  thereby  lift  the  degeneracy 
and  lower  its  energy  is  known  as  the  Jahn-Teller  effect. 


12 
In  the  more  general  case,  interaction  to  two  near-lying  states  via  a  suitable  nuclear 
motion  causes  the  surfaces  to  repel  each  other.  Looking  at  a  slice  of  the  two  surfaces 
along  the  symmetry-breaking  motion  that  couples  the  two  states,  the  upper  surface 
becomes  more  curved  as  the  interaction  increases,  while  the  lower  state  is  flattened. 
If  the  interaction  is  strong  enough,  the  lower  state  may  display  a  double  minimum  shape, 
with  the  original  high  symmetry  point  a  transition  state  between  the  two  minima.  This 
is  usually  referred  to  as  the  pseudo-Jahn-Teller  effect  (PJTE)  [7].  The  magnitude  of  the 
PJT  effect  depends  both  on  the  separation  of  the  two  states  and  the  magnitude  of  the 
linear  vibronic  constant. 

The  exclusion  of  linear  molecules  from  the  JT  effect  is  also  on  the  basis  of  symmetry 
arguments.  The  wavefunction  of  a  linear  molecule  is  either  even  or  odd  with  respect  to 
reflections  containing  the  axis  of  the  molecule  or  perpendicular  to  it,  so  of  course  the 
product  of  such  functions  will  always  be  even.  The  nuclear  motions  that  bend  the 
molecule  are  odd  in  this  respect,  and  so  can  never  be  contained  in  the  direct  product 
of  the  wavefunction  representations  —  the  linear  vibronic  constant  Fp  is  zero  for 

linear  molecules. 

For  linear  molecules,  and  other  cases  for  which  the  linear  vibronic  constant  is  zero 
or  negligibly  small,  the  quadratic  term  of  Eq.  2.14  comes  into  play.  For  linear  molecules 
in  degenerate  states,  the  quadratic  vibronic  coupling  leads  to  a  condition  similar  to  the  JT 
effect,  called  the  Renner-Teller  effect  (RTE)  [8,9].  Lifting  of  the  degeneracy  can  be  easily 
visualized  in  this  case.  Consider  a  linear  molecule  oriented  along  the  z  axis,  in  a  2U  state, 
so  that  one  component  of  the  tt  orbital,  say  ttx  contains  an  electron.   Bending  motions 
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in  the  orthogonal  directions  x  and  y  would  normally  be  equivalent  in  a  linear  molecule, 
but  if  one  component  of  the  it  is  occupied,  the  electronic  state  of  the  bent  structure  will 
be  different  depending  on  whether  the  electron  is  now  in  an  orbital  in  the  plane  of  the 
molecule  or  perpendicular  to  it.  Unlike  the  JT  case,  the  Renner-Teller  effect  does  not 
guarantee  that  lifting  the  degeneracy  by  bending  will  result  in  a  lower  energy.  The  two 
states  will  not  be  degenerate  (except  the  linear  geometry),  but  neither,  one,  or  both  may 
still  be  minima  at  the  linear  geometry,  or  the  linear  geometry  may  be  a  transition  state 
between  two  equivalent  bent  minima  depending  on  the  magnitude  of  the  interaction. 

In  the  more  general  case,  with  states  belonging  to  different  irreps  interacting,  there  is 
little  to  be  gained  from  a  symmetry  analysis  of  the  quadratic  vibronic  constant.  There  are 
usually  many  combinations  of  vibrational  motions  that  can  couple  two  electronic  states 
for  either  degenerate  or  nondegenerate  point  groups,  however  these  effects  are  usually 
relatively  small. 

Clearly  vibronic  coupling  effects  can  cause  major  problems  for  computational  meth- 
ods based  on  the  Born-Oppenheimer  approximation,  but  if  we  are  looking  for  minima 
on  potential  surfaces,  it  is  clear  that  such  points  are  usually  not  of  interest  —  when  the 
interactions  are  strong  enough,  they  cause  the  symmetry  of  the  molecule  to  be  broken, 
thus  moving  it  away  from  the  degenerate  or  pseudo-degenerate  point.  As  a  result,  it 
is  often  enough  to  recognize  points  on  the  PES  at  which  strong  vibronic  interactions 
occur  and  then  locate  and  follow  the  symmetry  breaking  to  lower  energies.  Of  course 
this  approach  will  not  always  work,  and  even  calculations  near  a  (pseudo-)  degenerate 
point  can  be  quite  taxing  for  the  method,  but  it  has  been  used  in  a  number  of  studies 
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(see  Bersuker's  bibliography  [10]  for  examples).  In  1984,  Lee  et  al.  [11]  demonstrated 
that  analytically  evaluated  hessians  could  be  used  to  compute  bending  frequencies  for 
Renner-Teller  molecules.  In  fact  if  due  care  is  taken  to  stay  on  the  correct  PES  during 
the  calculation,  hessians  evaluated  from  finite  differences  of  gradients,  or  even  energies, 
can  also  give  RT  bending  frequencies. 

2.6  Basis  Set  Expansion 

Several  approaches  may  be  taken  to  the  solution  of  the  electronic  Schrodinger 
equation,  but  the  most  widely  used  and  most  generally  applicable  is  expansion  in  a 
basis  of  Gaussian  functions.  These  functions  are  usually  (but  not  necessarily)  chosen 
to  mimic  the  atomic  orbital  (AO)  description  used  in  chemistry  and  are  referred  to  by 
that  name.  A  complete  (infinite)  basis  would  span  all  of  space  and  thus  allow  an  exact 
description  of  the  wavefunction.  In  practice,  however,  computational  resources  place 
limits  on  the  size  of  basis  that  may  be  employed,  and  it  is  necessary  to  compromise 
between  the  cost  of  the  calculation  and  the  accuracy  required. 

It  is  now  possible  to  evaluate  matrix  elements  of  the  Hamiltonian  operator  in  terms 
of  the  AO  basis  set.  For  a  set  of  functions  {\^},  we  have,  for  example,  the  one-  and 
two-electron  integrals  (or  matrix  elements)  [12] 


{n\h\v)  =  hpV  =    /  xl(xi)hxt(*i)dxi 

(fiu\Xa)  =   /  j£(Xi)x£(xa)  —  \A(x1)\(T(x3)dx1(ixa 
J  r12 


(2.16) 


15 

where 

.=  (*''     °r  (2.17) 

is  a  combination  of  the  position  and  spin  of  electron  i,  and  r\o  =  In  —  ro\.  The  left-hand 
sides  of  Eq.  2.16  are  shown  the  bra-ket  notation,  originally  due  to  Dirac  [13],  which 
is  in  common  use  and  for  the  one-electron  integral,  a  matrix  element  notation  is  also 
shown.  Both  the  bra-ket  and  matrix  element  notations  are  used  in  the  literature,  and  both 
will  be  used  in  this  work. 

2.7  One-electron  Approximations 

Given  a  basis  of  atomic  orbitals,  which  are  (very  rough)  approximations  to 
the  eigenstates  of  the  atom,  the  first  step  in  an  electronic  structure  calculation 
is  usually  to  find  a  set  of  one-electron  functions  for  the  molecule  as  a  whole, 
called  molecular  orbitals  (MOs).  These  orbitals  correspond  to  the  qualitative  ideas 
about  molecular  orbitals  often  used  by  chemists. The  most  common  prescription  for 
calculating  MOs  is  the  Hartree-Fock  (HF)  Self-Consistent  Field  (SCF)  procedure 
[12]. 

The  result  of  an  SCF  calculation  is  a  set  of  MOs  {<f)p},  which  are  represented  as 
linear  combinations  of  the  AO  basis  functions, 


=  Ec«av  <2-18) 


Jp 


with  the  matrix  C  determined  by  the  SCF  procedure.     Of  these  MOs,  some  will 
be  occupied  by  electrons   and  some  will  be   unoccupied,   or  virtual  orbitals.      The 
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antisymmetrized  product  of  the  occupied  orbitals  is  the  Slater  determinant  for  that 
state. 

Depending  on  the  particular  formulation  (primarily  the  treatment  of  electrons  with 
different  spins)  used  in  the  HF  procedure,  it  may  be  referred  to  as  RHF,  UHF,  or 
ROHF,  among  other  possibilities.  The  most  straightforward  is  restricted  Hartree-Fock 
(RHF),  applicable  only  to  closed  shell  molecules  (where  all  electrons  are  paired),  where 
the  alpha  and  beta  spin  orbitals  in  each  pair  are  restricted  to  have  the  same  spa- 
tial function.  The  restricted  open-shell  Hartree-Fock  (ROHF)  formalism  extends  the 
RHF  idea  to  open  shell  systems  by  requiring  that  all  paired  electrons  have  share  the 
same  spatial  function,  but  allowing  additional  orbitals  which  may  be  occupied  by  sin- 
gle unpaired  electrons.  In  the  unrestricted  Hartree-Fock  (UHF)  method  (also  known 
as  DODS  —  different  orbitals  for  different  spins),  the  alpha  and  beta  spin  electrons 
have  distinct  spatial  parts,  allowing  the  possibility  that  they  might  be  quite  differ- 
ent. 

Only  when  pairs  of  electrons  have  the  same  spatial  component  (RHF,  ROHF)  is  it 
guaranteed  that  the  HF  wavefunction  is  an  eigenfunction  of  the  S2  operator,  satisfying 
S2\HF)  =  s(s  +  l)\HF).  UHF  wavefunctions  may  suffer  from  spin  contamination, 
in  which  the  wavefunction  is  an  admixture  of  states  of  different  S2  values.  Because 
the  true  wavefunction  is  always  an  eigenfunction  of  Sa,  it  is  thus  possible  that  a  UHF 
wavefunction  may  be  to  some  extent  unphysical.  Nevertheless,  the  UHF  method  has 
proven  quite  useful,  particularly  as  the  basis  for  formulating  general  (i.e.  applicable  to 
both  closed  and  open  shell  systems)  correlated  methods. 
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2.8  The  Correlation  Problem 

Although  the  SCF  procedure  provides  a  very  useful  qualitative  description 
of  molecules,  it  is  generally  inadequate  for  quantitative  uses  requiring  high  ac- 
curacy. The  method  considers  each  electron  the  average  field  of  all  others, 
which  ignores  the  fact  that  the  motion  of  each  electron  is  instantaneously  cor- 
related with  all  others  because  of  the  Pauli  Principle.  For  the  higher  ac- 
curacy required  of  most  interesting  applications,  it  is  necessary  to  go  be- 
yond the  one-electron  approximation  and  treat  con-elation  effects  in  the  sys- 
tem. 

The  correlation  problem  is  usually  formulated  in  terms  of  the  interaction  between 
different  configurations  (occupations)  of  a  set  of  one-electron  functions.  The  molecu- 
lar orbitals  are  usually  used  for  this  purpose  because  the  idea  of  occupied  and  unoc- 
cupied orbitals  is  well  defined  and  the  derivation  of  equations  for  correlated  methods 
are  thus  greatly  simplified.  The  SCF  determinant  serves  as  a  reference  for  the  defi- 
nition of  the  excited  state  configurations  in  which  the  correlated  wavefunction  is  ex- 
panded. 

2.9  Occupation  Number  Representation 

The  occupation  number,  or  second  quantized,  formalism  [14]  is  a  convenient 
representation  for  the  development  of  correlated  methods  because  it  represents  in  a 
compact  way  the  changes  in  occupation  of  the  various  determinants.     In  this  de- 
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scription  and  throughout  this  dissertation,  the  convention  will  be  used  that  subscripts 
i,  j,  .  .  .  ,   n  refer  to  orbitals  which  are  occupied  in  the  reference  determinant,   a, 
b,  .  .  .  ,   f  refer  to  virtual  orbitals,    and   p,    q,  .  .  .  ,    u  refer  to  orbitals  of  either 
type. 

Creation  and  annihilation  operators  are  defined  that  act  to  create  or  destroy  an  electron 
in  a  particular  one-electron  function  (MO).  Thus  the  creation  operator  aP  would  create 
an  electron  in  orbital  p,  and  the  annihilation  operator  aq  would  remove  an  electron  from 
orbital  q.  The  Pauli  Principle  is  satisfied  by  the  fact  that  apap  =  aqaq  =  0,  which  is  a 
result  of  the  anticommutation  relations  among  the  creation  and  annihilation  operators. 

Strings  of  creation  and  annihilation  operators  act  on  the  reference  determinant,  |0), 
to  produce  excited  state  determinants.  For  example,  an  electron  can  be  promoted  from 
an  occupied  orbital  i  to  a  virtual  orbital  a, 

alat\0)  =  \*f)  =  |J)  (2.19) 

where  the  resulting  excited  state  has  been  shown  in  two  common  notations.  Similarly, 
higher  excitations  may  be  obtained 

aiakalajalai\0)  =  |*$)  =  |$)  (2.20) 


Operators  such  as  the  Hamiltonian  are  represented  in  the  second  quantized  formalism 
by  the  matrix  element  of  the  operator  in  the  basis  of  one-electron  functions  multiplied 
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by  a  string  of  creation  and  annihilation  operators: 


H  =  ^  hpga^aq  +  -   ^   {pq\\rs)a\)a\asar  (2.21) 

p,q  p,q,r,s 

where  we  have  used  the  antisymmetrized  two-electron  integral,  (pq\\rs)  =  (pq\rs)  — 


{pq\sr). 

Strings  of  creation  and  annihilation  operators  may  be  evaluated  using  Wick's  theorem 
[15]  and  the  concept  of  a  normal  order  of  operators.  Normal  ordered  operator  strings  are 
those  which  have  all  a]  and  all  aa  to  the  right  of  the  string,  which  can  be  accomplished 
using  the  anticommutation  relations  for  second  quantized  operators.  This  is  done  to 
exploit  the  fact  that  g||0)  =  aa\0)  =  0.  The  normal  order  of  an  operator  string  A  is 
denoted  N[A],  and  has  the  property 

(0|iV[A]|0)  =0.  (2.22) 

This  result  is  obvious  if  all  a]  and  aa  operators  are  brought  to  the  right,  where  they  act 
on  the  reference  giving  zero  immediately.  If  there  are  no  a\  or  aa,  then  the  operators 
act  on  the  reference  to  create  an  excited  state  on  the  left  hand  side,  which  is  orthogonal 
to  the  reference  function  on  the  right, 

.$?/'::)  =0.  (2.23) 

Note  too  that  in  the  absence  of  a\  or  aa,  A=N[A]. 

Wick's  theorem  tells  us  that  the  contraction  of  a  pair  of  operators  is  defined  as 


a^ag  -  N 


apal  =  N 


a\aq 


apa\ 


+  N 


N 


a\aq 


apa\ 


,  or 


(2.24) 
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with  the  following  definitions  according  to  whether  p  and  q  are  occupied  or  virtual  orbitals 
(they  must  be  the  same  type  since  Sai  =  0) 


N 


N 


N 


N 


a]  a. 


a,a 
1   3 


a\ab 


aaab 


Sij 


(2.25) 


=  0 


?ab 


Longer  strings  of  operators  are  evaluated  by  taking  all  possible  single  contractions,  then 
all  double  contractions  (two  simultaneous  contractions),  up  to  full  contractions,  with  the 
sign  in  each  case  determined  by  the  permutations  necessary  to  make  pairs  of  contracted 
operators  adjacent.  As  an  example,  consider  the  Hamiltonian,  Eq.  2.21.  The  one-electron 
term  is  cast  in  normal  ordered  form  by  recognizing  that 


a^cig  =  N 


N 


a\aq 


a]vaq 


+  N 


a\aq 


(2.26) 


+  S 


v 


so  that 


]T  hpqa]vaq  =  ]T  hpgN  a\aq\  +  J^  hi 


p.? 


p,y 


(2.27) 


The  two-electron  term  has  a  string  of  four  operators,  so  all  possible  single  and  double 
contractions  must  be  considered.  There  are  four  single  contractions  all  of  which  can  be 
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manipulated  into  the  same  form,  and  similarly  for  the  two  double  contractions, 


jY]  (pq\\rs)ala\a3ar  =  -   V   {pq\\rs)N 


p,q,r,s 


a^a]qasar 


1 


(2.28) 


+  J3(p»lk*)JV  alai  +oX)(i.2lK?) 


P,Q,i 


»J 


When  evaluating  a  product  of  several  normal  ordered  operators,  the  same  rules  hold,  but 
all  contractions  which  are  internal  to  any  given  operator  have  already  been  accounted 
for,  so  only  those  contractions  connecting  different  normal  ordered  operators  need  be 
considered. 

When  developing  correlated  methods,  it  is  convenient  to  employ  the  normal  ordered 
product  form,  Oyv,  of  a  second  quantized  operator,  0,  defined  by 


0N  =  O  -  (0|O|0>. 


(2.29) 


This  form  has  the  property. 


<0|Ojv|0>  =  (0|O  -  (0|O|0)|0)  =  (0|O|0>  -  (0|O|0)  =  0 


(2.30) 


which  allows  one  to  focus  on  the  correlation  correction  to  the  result  rather  than  carrying 
the  reference  expectation  value  (a  constant)  through  all  of  the  equations.  For  example, 
consider  the  Hamiltonian.  The  normal  ordered  operator  is  the  sum  of  Eqs.  2.27  and 
2.28,  or 


N[H]  =  J2hPiN[4^]+1EPiq^Pi^ 


N 


a\aq 


1 


+  tY]        (pq\\rs)N 

a    <     'Ti.a.r.s 


4  -^^p.g,r,s 


o-la\asar 


hj 


(2.31) 
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and  the  expectation  value  is 


(0|tf|0)  =  (0 


Y^  hpqa]paq  +  -  J2p,q>r,s  (Pl\\rs)alala*a' 


=  E>i»+?E< 


(2.32) 


Ull*J 


»j 


The  normal  ordered  product  form  of  the  Hamiltonian  is  thus 
HN  =  H-(0\H\0) 

=  Y/hpqN\alaq]+Toi(P^\\^)N\alag]+jTors(P({\\rs)N 


p,q 


'p,q,t 


<rpaqasar 


=  E/«Jvk«t]  +  iE..r.(MiiM>JV 


P. 9 


ajajasar 


(2.33) 


with  the  Fock  operator, 


In  =  E  fp*N  [4a*]  -  E  ^  [°J°«]  +  E  <P»ll9»)^ 


aja9 


(2.34) 


p,q 


p,? 


p,?,» 


used  in  the  final  form.  The  normal  ordered  product  Hamiltonian  is  also  frequently  written 
in  terms  of  its  one-  and  two-electron  components, 


HN  =  fN  +  WN. 


(2.35) 


CHAPTER  3 
COUPLED  CLUSTER  METHODS 


3.1   Introduction 

Two  of  the  most  common  approaches  currently  in  use  for  obtaining  correlated  energies 
and  properties  are  configuration  interaction  (CI)  and  coupled  cluster  (CC)  methods. 
Both  methods  use  an  expansion  space  of  excited  Slater  determinants  to  represent  the 
wavefunction;  the  difference  between  them  originates  with  the  particular  form  chosen  for 
the  expansion.   In  the  CI  case,  a  linear  expansion  is  used, 

\q>cl)  =  (l  +  C)|0)  =  (i  +  ft  +  C2  +  •  •  -)|0)  (3.1) 

where  C,  generates  all  possible  i-fold  excitations  from  the  reference  determinant  weighted 
by  coefficients  giving  the  importance  of  each  determinant  in  the  total  CI  wavefunction. 
In  most  cases,  it  is  necessary  for  computational  reasons,  to  truncate  the  expansion  space 
somewhere  short  of  including  all  possible  excitations  (up  to  rc-fold  for  an  n  electron 
system).  This  results  in  a  variationally  determined  wavefunction  (and  energy),  but 
sacrifices  the  size  extensivity  property. 

Size  extensivity  [2,3]  is  the  requirement  that  the  energy  of  a  system  scale  properly 
with  the  number  of  electrons.  It  is  most  easily  understood  by  considering  an  assembly  of 
A^  identical  non-interacting  molecules.  If  each  molecule,  p,  has  an  energy  and  reference 
function  satisfying 

//|0(p))  =  £#°|0W)  (3.2) 
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(the  one-electron  functions  of  |0(?'')  are  assumed  to  be  localized  on  that  molecule),  then 
we  can  write  the  reference  energy  and  wavefunction  of  the  entire  assembly  as 

E(total)  =  J2EiP) 

P=l 

(3.3) 

N 

|0(*>*«0)  =  Yl  |0(p)) 
p=l 

Since  the  molecules  are  non-interacting,  we  expect  the  correlated  energy  and  wavefunc- 
tions  to  be  representable  as  similar  sums  and  products,  respectively.  When  a  truncated  CI 
expansion  is  applied  to  this  system,  however,  this  is  not  so.  For  an  individual  molecule, 
p,  the  coefficients  C\ |  ,  C^  ,  •  •  •  Cm  .  with  m<n,  are  determined  by  the  Schrodinger 
equation, 

H  (i  +  c[p)  +  Cf  +  ■  ■  ■)  |0(p))  =  41  (l  +  C[p)  +  C[p)  +  ■■■)  |0<p»)  (3.4) 

and  give  a  product  wavefunction  which  includes  excitations  up  to  Nm  due  to  products 
of  excitations  on  different  molecules.  If  the  assembly  were  treated  as  a  whole,  however, 
only  up  to  7?i-fold  excitations  would  appear  in  the  wavefunction.  Thus,  they  can  never  be 
equivalent.  The  energy  of  the  assembly  suffers  similar  problems.  For  example,  a  lattice 
of  N  non-interacting  hydrogen  molecules  treated  with  CI  restricted  to  double  excitations 
(CID)  has  an  energy  which  scales  as  \/N  instead  of  N,  resulting  in  an  error  of  12.3%  for 
10  molecules  and  21.1%  for  20  [16].  Calculations  by  Ahlrichs  and  coworkers  at  the  CID 
level  show  that  the  difference  between  supermolecule  (the  two  molecules  separated  by  a 
very  large  distance  calculated  together)  and  the  sum  of  individual  molecular  energies  are 
~7  kcal  mol"1  for  BH3  +  BH3  [17],  and  -15  kcal  mol"1  for  CH3F  +  P  [18]. 
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The  well  known  property  of  exponentials,  eAeB  =  eA+B  (where  [A,B]=0),  suggests 

an  alternative  solution  to  the  problem.  If  the  individual  molecule  wavefunctions  are  of 

the  form 

|*W)  =er<P)|0(p)}  (3.5) 

with  T"  being  an  expansion  operator  analogous  to  C,  in  Eq.  3.1,  and  satisfying 

He7<P)\oM)  =  E^e'I<P)\0^)  (3.6) 

then  the  product  wavefunction  is 

N 
\y(total)j  =  TTeT<p)|0(7))^ 


p=l 


r^    n  (3.7) 

n  i°w> 


p 


=  er<'-"|0(/o/a/)) 
which  is  identical  to  the  total  wavefunction  resulting  from  a  calculation  of  the  entire 
assembly.  This  is  the  form  of  the  wavefunction  employed  in  coupled  cluster  theory. 

This  simple  physical  picture  demonstrates  the  need  for  the  exponential  wave  operator, 
but  it  is  the  Brueckner-Goldstone  linked  diagram  theorem  [19-23],  which  CC  methods 
satisfy  [24],  that  guarantees  size  extensivity  when  this  form  is  used.  Thus  when  the  system 
is  composed  of  interacting  fragments,  the  simple  physical  picture  becomes  clouded,  but 
the  size  extensivity  property  is  still  guaranteed. 

Having  motivated  the  use  of  coupled  cluster  methods  and  the  exponential  wave 
operator,  in  the  remainder  of  this  chapter,  we  will  derive  the  working  equations  of 
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coupled  cluster  theory  and  other  results  which  will  be  needed  for  the  Fock-space  CC 
methods  discussed  in  the  next  chapter. 

3.2  Basic  Equations 

In  order  to  expand  the  coupled  cluster  wavefunction  in  the  space  of  Slater  deter- 
minants, we  define  a  set  of  cluster  operators  which  are  capable  of  representing  the 
contributions  to  the  wavefunction  due  to  determinants  of  particular  levels  of  excitation, 

T\  =  ^  t*a\a.i 
i,a 

T2  =  ^2  t^a\a\aiaj 

(3.8) 
T3=   Y^  *SaMaca«'aia* 


>>J>k 

a>b>c 


where  the  cluster  amplitudes,  ^*£ ..'  are  the  scalar  coefficients  for  each  excited  Slater 
determinant  in  the  wavefunction.  The  cluster  operators  for  each  level  of  excitation  are 
collected  together  into  a  single  operator, 

T  =  T!  +  T2  +  Tz  +  ■  •  •  (3.9) 

for  convenience. 

To  obtain  the  multiplicatively  separable  coupled  cluster  wavefunction,  the  wave 
operator  (which  produces  the  correlated  wavefunction  from  the  reference)  is  exponential 
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in  the  cluster  operators, 

|*cc)=ftcc|0)  =  er|0).  (3.10) 

This  wavefunction  must  satisfy  the  Schrodinger  equation 

HeT\0)  =  EcceT\0)  (3.11) 

where  E„  is  the  coupled  cluster  energy  and  H  is  the  second  quantized  electrostatic 
Hamiltonian.  In  order  to  simplify  later  discussion,  we  change  to  the  normal  product  form 
of  H  by  subtracting  Eq  =  (0\H\0)  from  both  sides  of  Eq.   3.11, 

(H  -  (0\H\0))eT\0)  =  (Ecc  -  Eo)eT\0) 

(3.12) 
HNeT\0)  =  EcorreT\0) 
where  the  correlation  contribution  to  the  energy  has  been  designated  E„„.   Finally,  we 

multiply  both  sides  of  Eq.    3.12  by  e~T  to  obtain 

e-THNeT\0)  =  ECOTT\0)  (3.13) 

Notice  that  the  operator  e~TH^eT  acts  as  an  effective  Hamiltonian  operator,  acting  on 
the  reference  state  to  give  the  con-elation  energy;  this  quantity  is  often  referred  to  as  H^ 
in  post-CC  methods  (i.e.  those  for  which  the  derivation  relies  on  having  first  solved  for 
a  CC  energy  and  wavefunction). 

The  effective  Hamiltonian  can  be  simplified  by  the  use  of  a  Baker-Campbell- 
Hausdorff  expansion 

e~THeT  =  H+[H,T]  +  h[H,T],T]  +  ---  (3. 14) 
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The  commutators  of  strings  of  second  quantized  operators  may  be  evaluated  using  the 
contraction  theorem  of  Harris,  Jeziorski  and  Monkhorst  [25],  with  the  result  that  the 
expansion  terminates  after  the  quadruply-nested  commutator  term  because  H  (Eq.  2.21) 
contains  only  four  operators.   This  result  is  often  denoted 


e~THeT 


=  (//eT)f  (3.15) 


where  the  subscript  c  indicates  that  only  so-called  connected  terms  remain.  The  ter- 
minology arises  from  the  diagrammatic  formulation  of  the  equations  often  used  in  CC 
theory:  surviving  terms  are  those  in  which  their  diagrammatic  representation  forms  a 
connected  figure,  from  which  no  component  may  be  removed  except  by  breaking  lines. 
The  algebraic  analog  is  that  the  term  cannot  be  separated  into  a  product  of  distinct  sum- 
mations. Disconnected  terms  are  those  which  can  be  separated  into  two  components 
without  breaking  a  line,  or  algebraically  the  term  may  be  cast  as  the  product  of  two 
distinct  summations.  Use  of  Eq.  3.15  greatly  simplifies  reduction  of  the  operator  form 
of  the  coupled  cluster  equations  to  a  spin  orbital  form  suitable  for  numerical  work. 

The  equations  which  determine  the  CC  energy  and  amplitudes  are  obtained  from  Eq. 
3.13  by  projecting  from  the  left  with  the  reference  determinant  and  each  class  of  excited 
determinants  (single,  double,  triple,  etc.).  These  projections  provide  equations  which 
determine  the  energy,  T„  To,  T3,  etc.,  respectively. 

{0\^HNeT>j  \0)  =  Ecorr  (3.16) 
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a 
i 

ab 
ij 

abc 
ijk 


K), 

M, 


0)  =  o 


0  )  =0 


(3.17) 


3.3  Reduction  to  Spin  Orbital  Equations 

The  operator  expression  for  the  energy,  Eq.  3.16,  can  be  expanded  to 

o -     ~     '-'    -     --     1 


0)     (3.18) 


(In  +  WN)ll  +  Tx  +  T2  +  -T{  +  T3  +  TxT>  +  -T{  +  • 

by  replacing  Hx  with  its  one-  and  two-body  components  (see  Eq.  2.21)  and  expanding 
the  exponential.  Since  the  Hamiltonian  is  in  normal  product  form,  the  lead  term  in  the 
exponential  expansion  is  (0|//#|0)  =  (0|(//  —  (0|//|0))|0)  =  0.  The  remaining  cluster 
operators  act  on  the  reference  to  give  an  excited  determinant  and  a  product  of  amplitudes. 
As  a  one-body  operator,  /#  can  create  at  most  a  singly  excited  state  acting  to  the  left  on 
the  bra  reference  determinant,  consequently  nonzero  contributions  will  result  only  from  a 
singly  excited  state  on  the  right.  Similarly,  the  two-body  W^  operator  requires  a  doubly 
excited  state  on  the  right.   Thus,  the  equation  simplifies  to 


(ra  +  itf)|o) 


fNTi  +  WN[T2  +  ^Tl 


(3.19) 
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It  should  be  observed  that  the  CC  energy  includes  contributions  due  only  to  T, 
and  T2  regardless  of  higher  level  cluster  operators  in  T;  higher  clusters  contribute  only 
indirectly  via  their  effects  on  Tt  and  T2. 

Reduction  of  the  amplitude  equations  (Eq.   3.17)  to  a  more  explicit  operator  form 
proceeds  in  a  similar  fashion.  For  example,  the  7,  equation  can  be  expanded  to 


0  = 


(  "  Un  +  wN)  (i  +  r,  +  t2  +  Ii?  +  n  +  T1T2  +  It?  +  ■  ■  •)  o\     (3.20) 


A  nontrivial  equation  is  obtained  when  the  combined  action  of  Hs  and  the  cluster 
expansion  produces  a  single  excitation  |f)  on  the  right.  The  one-body  operator  fN 
is  capable  of  changing  the  excitation  level  by  -1,  0,  or  +1,  and  WK  by  -2,  -1,  0,  +1,  or 
+2.  The  fN  +1  term  and  the  WN  +2  term  cannot  form  connected  contractions  with  T,  so 
they  can  only  contribute  to  the  T,  and  T2  equations,  respectively,  via  the  "1+"  term  in 
the  exponential  expansion.  Because  of  the  nature  of  the  Hamiltonian,  mentioned  above, 
it  should  be  evident  that  the  equation  for  a  particular  cluster  amplitude  Tm  will  involve 
contractions  of  fN  with  products  of  cluster  operators  resulting  in  m  and  rn+1  excitations, 
and  WN  with  cluster  operator  products  yielding  m-1  to  m+2  excitations.  Equations  for 
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the  first  three  cluster  operators  are  thus 


fN[l  +  Ti  +  T2  +  yTl 


T1  +  T2  +  ^Tl  +  Ti  +  T\T2  +  ^T[ 


0  =/°.b  fN{T2  +  T3  +  TiTj}  +  WN  j  1  +  Ti  +  T2  +  ^T'i  +  Tz  +  T^ 


0  = 


+^T?  +  T4  +  TlT,  +  ^Tl+  ±T?T2  +  ^T? 

abc 
ijk 


fNiTA  +  T4  +  TlT1+^ 


T2  +  T3  +  TXT2  +  T4  +  T.Tz  +  ir22  +  ^T{T2 


+T5  +  T,T4  +  ir27s  +  7^73  +  ^T,Ti  +  ^T?T2 


(3.21) 


where  the  operator  determined  by  the  equation  has  been  underlined. 

Conversion  of  the  operator  equations  into  the  spin  orbital  equations  useful  for 
numerical  work  is  done  by  considering  in  turn  each  of  the  operator  products  in  the 
foregoing  equations  and  using  Wick's  Theorem  [15]  to  evaluate  the  possible  contractions 
among  the  strings  of  second  quantized  operators  contained  in  the  Hamiltonian  and  cluster 
operators.  This  analysis  can  be  carried  out  using  a  purely  algebraic  approach  or  a 
diagrammatic  approach  [26],  but  rather  than  adding  greatly  to  the  length  and  complexity 
of  this  overview  we  prefer  to  omit  the  detailed  procedure  here,  and  simply  report  the 
result. 
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Consider,  for  example,  the  first  term  of  the  CC  energy  as  given  by  Eq.  3.19: 

(o|/jvri|o>=  E  fP^(o\{4ai}{ala'}\°] 


p,q,i,a 


=  E  /w*?M««<o|o) 

p,q,i,a 

=  £/«*■ 


(3.22) 


i, a 


Similarly  for  the  second  term. 


T2  +  \Tl 


0     = 


p,q,T,s,i,j,a,b 


;s,i,j,a,b 

0  {aJa{a*0r}({aWa*°i}  +  {aiai}{aIai}) 
i       E       (Mil")  (#  +  |tf^Mi>4*M0|0) 


■p,q,r,s,i,j,a,b 


i,j,a,b 


(3.23) 


which  gives  the  spin  orbital  equations  for  the  coupled  cluster  energy: 


ECOrr  =  E  /»*?  +  J    E    (Vllfl6)  (*«   +  ^*i)  (124) 

i,a  i,j,a,b 


Spin  orbital  forms  of  the  amplitude  equations  are  obtained  in  the  same  way.  Complete 
spin  orbital  equations  for  the  CCSDTQ  model,  the  most  complete  CC  model  implemented 
at  the  present  time,  may  be  found  in  the  recent  paper  of  Kucharski  and  Bartlett  [27]. 
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3.4  Numerical  Solution  of  the  Coupled  Cluster  Equations 

The  overall  result,  as  can  be  seen  in  Eq.  3.21  is  a  set  of  coupled  nonlinear  equations 
for  the  amplitudes  £",*£"•'  based  on  the  MO  integrals  fpq  and  (pq\\rs).  These  equations 
are  usually  solved  by  fixed-point  iteration  starting  with  a  guess  for  t^  obtained  from 
second-order  perturbation  theoiy,  and  it  has  been  found  that  acceleration  methods  based 
on  conjugate  direction  techniques  [28-32]  are  very  useful  in  obtaining  faster  convergence, 
and  thus  reducing  the  computational  effort.  Once  the  amplitudes  are  converged,  Eq.  3.24 
may  be  evaluated  to  obtain  the  CC  energy.  In  many  cases,  the  amplitudes  and  energy 
will  be  used,  along  with  the  MO  integrals,  as  inputs  to  further  calculations  to  determine 
the  CC  density  matrix,  gradients  of  the  energy,  or  other  results,  for  which  it  is  often 
useful  to  work  in  terms  of  the  CC  effective  Hamiltonian. 

3.5  The  Coupled  Cluster  Effective  Hamiltonian 

When  the  coupled  cluster  calculation  is  carried  out  as  the  first  step  in  a  Fock-space 
CC  calculation,  it  is  useful  to  form  the  elements  of 

HN  =  e~THNeT  =  (HNeT^  (3.25) 

explicitly  from  the  MO  integrals  and  CC  amplitudes  because  the  FSCC  equations  can 
be  expressed  in  a  compact  and  computationally  efficient  form  using  H^.  Note  that 
this  operator  is  non-Hermitian,  so  the  bra-ket  symmetry  of  the  MO  integrals  (i.e. 
(p§||rs)  =  {r±\\pq))  is  not  carried  over  to  the  elements  of  Hfir:(pq\\rs)  ^  (rs\\pq). 
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Because  //jv  contains  at  most  four  second  quantized  operators,  H^  will  include 
contractions  of  H^  with  as  many  as  four  T  clusters.  Like  H^,  the  effective  Hamiltonian 
is  described  as  having  one-  and  two-body  components.  In  fact,  while  Hm  terminates 
at  two-body  terms,  Hn  has  up  to  n+2-body  terms  in  an  n  electron  system.  Just  as  the 
Hamiltonian  is  often  written  in  terms  of  /#  and  W^,  the  components  of  the  effective 
Hamiltonian  are  usually  denoted  /#,  W^,  II I /v,... 

The  Hj\r  elements  required  for  the  Fock-space  coupled  cluster  methods  considered  in 
the  next  chapter  include  all  of  the  one-  and  two-body  terms  and  a  few  of  the  three-body 
terms.  They  are  obtained  in  the  same  manner  as  the  spin  orbital  equations  for  the  CC 
energy  and  amplitudes,  above,  and  the  results  are  shown  below. 

The  first  stage  is  to  isolate  the  components  of  the  expansion 


(HNeT)c  =  Ufa  +  WN)(l  +  li  +  T2  +  ^T{  +  T3  +  T,T2  +  ilj  +  •  •  ^j) 

(3.26) 
from  which  contractions  can  form.  At  the  operator  level,  the  equations  are  closely  related 
to  the  CC  amplitude  equations  (Eq.  3.21),  but  with  two  important  differences.  First,  H^ 
is  defined  based  on  converged  cluster  amplitudes,  so  Eq.  3.21  must  be  satisfied  prior  to 
forming  the  effective  Hamiltonian  elements.  Secondly,  like  the  Hamiltonian  itself,  the 
effective  Hamiltonian  can  have  -n,  .  .  .  ,  0,  .  .  .  ,  n  excitation  character  for  an  n-body 
component.  So  while  ("£'''  |///v|0 \  =  0  is  identical  to  the  CC  amplitude  equations  given 
above,  the  effective  Hamiltonian  itself  contains  numerous  terms  besides  those  required 
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to  satisfy  the  equation.  Thus  the  one-body  component  is 

fN=(fNll  +  Tl+T2  +  ^ 


+wN{  t1  +  t2  +  irj2  +  n  +  t,t2  +  irf 


c 


and  the  two-  and  three-body  components  of  H^,  are 

WN  =(fN{T2  +  T3  +  TXT2}  +  WNll  +  T,+T2  +  ITj2  +  T3  +  T{T2 

+^r13  +  t4  +  T&  +  ^r22  +  ^t?t2  +  ir,4 


and 


h  l<T2Tz  +  ^T{T3  +  ^TiT22  +  ^T?T2\  J 


2!  2!    1  2! 


(3.27) 


(3.28) 


+  wN  j  t2  +  t3  +  Ti  r2  +  r4  +  Ti  r3  +  ^r22  +  ^r^  (3.29) 


Reduction  to  spin  orbital  form  follows  the  same  procedure  used  for  the  CC  equations. 
In  general,  the  /z-body  component  of  the  effective  Hamiltonian  will  have  (n+1)2  different 
terms  according  to  the  distinct  labelings  of  the  indices  as  occupied  or  virtual  orbitals. 
The  one  corresponding  to  the  n-fold  excitation  will,  of  course,  be  zero  and  for  higher  n, 
an  increasing  number  of  other  terms  will  be  zero.  For  example  in  III  n  it  is  impossible 
to  form  terms  with  -3,  -2,  or  -1  excitation  character,  as  well  as  half  of  the  possible 
0  excitation  terms  (the  all  occupied  and  all  virtual  labelings).  Thus  while  the  number 
and  complexity  of  higher-body  effective  Hamiltonian  terms  does  increase,  it  is  not  as 
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explosive  as  that  seen  in  many-body  perturbation  theory  (MBPT)  [26].  The  spin  orbital 
equations  for  H^  are  give  in  Tables  3.1-3.3. 

Table  3.1.  One-body  coupled  cluster  effective  Hamiltonian  elements  in  spin  orbital  form. 
Excitation  level  of  operator  is  given  in  the  first  column.  See  text  for  explanation  of 
notation. 


te 


te 


te 


-1      fal=  fia  +  {im\\ae)ter 

0     k  =  fij  +  fiftf  -  {im\\je)tem  +  i(*m||/c)t£,  +  (im\\fe)tf^ 

0     Lb  =  fab  -  fmbC  +  (am\\be)tem  -  \{nm\\be)t^m  -  (nm\\be)tanter 

■"■'        Jia  =  Jai        Jmi^rn    '    Ju-e'i     '    Jmeiim        Jine'^m        \m(l\re)tjn 

-\(nm\\ie)tanem  -  (nm\\ie)tantem  +  |(am||/e}*£  +  {am\\fe)t{t'm 

+Iv*H|e/)C»  +  (mn\\ef)t^tfn  -  ±(m,n\\ef)tpaJn 
-k(mn\\ef)tltl[  -  (mn\\ef)t^ntfn 


Table  3.2.  Two-body  coupled  cluster  effective  Hamiltonian  elements  in  spin  orbital  form. 
Excitation  level  of  operator  is  given  in  the  first  column.  See  text  for  explanation  of 
notation. 


-2      (ij\\ab)  =      (ij\\ab) 
-1      (ij\\ka)=     (ij\\ka)  +  (ij\\ea)tl 
-1      {ai\\bc)  =      (ai\\bc)  -  (mi\\bc)tam 
0      (ij\\kl)  =      (ij\\kl)  +  {-fP{kl){ij\\ke)t\  +  \{ij\\ef)4 

+(-)PP(kl)(ij\\ef)tlt{ 
0      (ab\\cd)  =      (ab\\cd)  -  (-)pP{ab)(am\\cd)thm  +  \_(mn\\cd)t^n 

+(-)pP(ab)(mn\\cd)tamtbn 
0      (ia\\jb)  =     (ia\\jb)  -  {im\\jb)tam  +  (ia\\tb)t)  -  {im\\eb)tfm  -  (im\\eb)tp"m 
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Table  3.2.  (Continued)  Two-body  coupled  cluster  effective  Hamiltonian  elements  in 
spin  orbital  form.  Excitation  level  of  operator  is  given  in  the  first  column.  See  text  for 
explanation  of  notation. 

+  1      (taffjAr)  =     (ia\\jk)  +  fiet™  -  (im\\jk)tam  +  (-)pP(jk)(ia\\je)tl 
+(-)pP(jk)(im\\je)t%k  -  (-YP(jk)(im\\je)tttam 
+(-YP(jk)(ia\\ef}tp{  +  \{ia\\ef)t%  +  l{mi\\ef)t%k 
+<^||e/)^-i(im||/e)^C 

-(-)pP(jk)(im\\fe)tf£tl  -  (-)pP(ifc)(im||/e)^|^ 
+  1       (ab\\cr)  =      (ab\\ci)  -  fmctti  +  (ab\\ce)t\  -  (-)" P(ab)(am\\a)tbm 
+  ±(mn\\cr)t±  +  (-)pP(ab)(mn\\ci)tamtbn 
+(-)pP(ab){am\\ce)t%t  -  {-)p  P{ab)(am\\ce)t\thm 
-i("»n||cc>Cti  "  (mn||ec>^t«{  +  (-)*P(a&)(nm||ce)i£i?4 
+  ±(nm\\ce)t±n  +  (-)PP(a&)<nm||ce)i^4 
+2      (fl6||y)  =      (a6||»j)  -  (-)pP(ij)fmjt£  +  (-)pP(ab)fbet% 

+/«.*&  -  {-)PP(ij)fmettt)  ~  (-)pP(ab)fmet%tbm 
-(-rP(ab)(mb\\ij)t«m  +  (-)pP(ij)(ab\\ej)t<  +  \{mn\\i])t± 
+  \(ab\\ef)teJ  -  (-)pP(ij)(mb\\ie)t%3  +  (-)p P(ab)(mn\\tj)tamtbn 
+(-)pP(ij)(ab\\ef)titfj  -  (-)pP(ij)<m6||ie)C*S 
-\{-YP{iJ)[nm\\3e)ttm  +  \{-?P{ab)(bm\\fe)^ 
-(-)pP(ii)(mn||ei)^^.  +  (-fP(a6)(ma||e/)^if/ 
-(-)pP(*i)(a6)(nm||ic>*-<en;j ■  +  (-Wi)(a&)(am||/e)t{t*- 
+i(-)pP(u)(mn||ej)^C6n  -  |(-)PP(a&)(m&||e/)Oi/ 

+(-M»i)(a6)(mn||Cj)*JtJJl4-(-)'P(y)(a6)(m6||e/)*X*i 
+\(mn\\ef)tf°Jn  +  (mn\\ef)tfjn  -  i(-)"P(u)(mn||e/)^ 

4(-)pP(a6)(mn||e/)C*-5  +  (-Mv)(mn||e/)C*3 
+J<mn||c/>t?/tS*fi  -  K-)pP(a6)(mn||e/)if/^n 


38 

Table  3.2.  (Continued)  Two-body  coupled  cluster  effective  Hamiltonian  elements  in 
spin  orbital  form.  Excitation  level  of  operator  is  given  in  the  first  column.  See  text  for 
explanation  of  notation. 


-I(-)pP(ij)(mn||e/)t^-(-fP(ii)(a&){mrl||e/)tX^ 

+  ^-)PP(ij)(mn\\ef)t<tfJt±l  +  \{-f  P{ab){mn\\ef)tamthnt\j 

-(-)pP(ij)(nm\\fe)tfntp%3  -  (-)pP(ab)(nm\\fe)tfntamt^ 
+  (-fP(lJ)(mn\\ef)txMtbn 

The  summation  convention  (all  repeated  indices  are  summed)  has  been  used  to  save 
space  and  simplify  presentation  of  these  formulae.  The  permutation  symbols,  (  —  )PP(), 
are  used  in  the  same  fashion  as  Lee  et  al.  [33]  and  Kucharski  and  Bartlett  [27].  A 
permutation  of  two  indices,  e.g.  (  —  )pP(ij),  denotes  two  terms.  One  is  the  identity 
permutation,  with  signature  p=+l,  while  the  second  is  the  term  with  indices  i  and  j 
interchanged.  A  permutation  of  three  indices,  e.g.  (  —  )pP(ijk),  represents  six  terms  with 
signature  p=+l  for  cyclic  permutations  and  p=-l  for  acyclic  permutations.  The  symbol 
(  —  )pP{i/jk)  indicates  that  permutations  i  <->■  j  and  i  <->  k  should  be  taken,  but  not 
j  <->  k;  the  signature  is  the  same  as  for  a  binary  permutation.  Finally,  a  composite 
symbol,  such  as  (—)pP(i/jk,ab),  indicates  that  the  direct  product  of  the  two  sets  of 
permutations  (i/jk  and  ab)  should  be  taken. 

3.6  Truncation  of  the  Cluster  Operator 

It  is  worth  noting  that  we  have  made  no  approximations  in  deriving  the  equations 
given  above,  and  the  energy,  wavefunction,  and  effective  Hamiltonian  thus  determined 
would  be  exact  (for  the  chosen  atomic  orbital  basis).    Unfortunately,  the  number  of 
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Table  3.3.  Three-body  coupled  cluster  effective  Hamiltonian  elements  needed  for 
Fock-space  (1,0),  (0,1),  and  (1,1)  sectors  in  spin  orbital  form.  (The  full  IIIjv  includes 
an  element  at  excitation  levels  +2,  and  +3  in  addition  to  those  shown  here.)  Excitation 
level  of  operator  is  given  in  the  first  column.  See  text  for  explanation  of  notation. 

0     (aij\\klc)  =    (ij||cc)«Jf 

0      (abi\\jcd)  =   -(mi\\cd)tfm 
+  1      (ija\]klm)  =  (-)PP(ml/k)(ji\\ek)t^  +  ±(j*||e/)Ci 

+(-)PP(k/lm)(ji\\ef)t^tfk 
+  1       (abc\\dei)  =    -(-fP(b/ac)(m.b\\de)t™  +  \(rnn\\de)t^n 

M-)PP(b/ac)(mn\\dc)tZthn 

+  1      (iab\\jck)  =    -(mi\\ec}1%  +  (-)pP(jk)(mi\\ec)tfmtl 

+(-yP(ab){mi\\ec)tpbm  -  (-)pP(jk)(mi\\jc)tlbm 
+  {-)* '  P{ab){ia\\ec)tbk) 

+2      (iab\\jkl)  =    (-)PPUk/l)(im]\jk)t^t  +  (-)pP(j/kl:ab)(ia\\je)trt 

-(-)pP(A;/j7,afe)(mi||fce)0^-(-)pP(ifc/)(im||efc)^ 
+(-)pP(jA://,a6)(im||/e)^^-i(-)pP(jA://)(im||e/)ieJ^ 
+(-)pP(j/kl,ab)(ia\\ef)tfkb  +  /iet$ 
+(-)pP(j/kl)(im\\je)t%bkl  +  i(-?P{ab)(ia\\ef)t$ 
M-rPU/klHimllfe)^  -  \{-YP{ab)(im\\ef)tak){tbm 

+(-)'i'O'^a*)(«m||/e>^*i-l(-)'i'(;*/0<*"»lk/>*Sf*i 

+  ^m»l|e/)^>  +  (m»||e/)^/t 

amplitudes  to  be  determined  is  0(N!)  for  N  AO  basis  functions.  This  rapidly  becomes  too 
expensive  to  compute  using  current  and  forseeable  computer  facilities,  so  it  is  necessary 
to  make  approximations.  The  obvious  choice  is  to  limit  the  cluster  operators  included 
in  T  to  lower  levels  of  excitation  because  the  number  of  amplitudes  grows  by  a  factor 
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of  n(N-n)  for  each  additional  excitation.   Because  of  the  exponential  nature  of  the  CC 
wave  operator,  some  higher  excitation  effects  are  incorporated  by  products  of  lower-level 
excitations  (i.e.  T,  T2  is  a  triple  excitation).  In  fact  in  most  cases,  such  products  dominate 
the  wavefunction  at  a  given  level  of  excitation  [34]. 

The  single  excitations  alone  cannot  have  any  contribution  to  the  energy,  so  the  first 
reasonable  approximation  is  T=T2,  known  as  the  CCD  (coupled  cluster-doubles)  model. 
With  the  doubles  present,  the  singles  can  also  contribute,  so  T=T1  +  T2,  the  CCSD 
model,  would  be  next  on  the  hierarchy.  Increasingly  sophisticated  methods  are  obtained 
by  adding  the  next  level  cluster  operator:  CCSDT,  CCSDTQ, 

To  obtain  approximate  versions  of  the  equations  given  in  this  chapter,  simply  consider 
all  cluster  operators  not  present  in  the  desired  model  to  be  zero  and  eliminate  them  from 
the  equation.  For  example,  CCSD  equations  would  have  T3=T4  =  T5=0,  and  all  terms 
involving  these  operators  would  be  eliminated.  Observe,  however,  that  products  of  Tt 
and  T2  giving  overall  triple,  quadruple,  and  higher  excitations  would  still  appear.  This  is 
a  consequence  of  the  exponential  expansion  and  is  the  reason  why  a  given  truncation  of 
coupled  cluster  theory  recovers  more  of  the  correlation  energy  than  the  same  truncation 
of  operator  manifold  in  the  linear  CI  expansion  (Eq.  3.1). 

In  practice,  the  CCSD  model  has  proven  widely  applicable,  and  the  more  expensive 
CCSDT  model,  and  approximations  which  involve  perturbative  inclusion  of  the  triple 
excitation  effects,  provide  higher  accuracy  and  have  been  shown  capable  of  handling 
some  problems  that  would  otherwise  require  multireference  techniques  [34].  The  CCSDT 
model,  or  perturbative  approximations,  are  currently  the  practical  limit  for  the  majority 
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of  interesting  systems.  Further,  the  computer  code  used  in  this  work  does  not  presently 
include  the  full  CCSDT  model,  so  calculation  of  Hn  is  limited  to  the  CCSD  level. 


CHAPTER  4 
FOCK  SPACE  COUPLED  CLUSTER  THEORY 


4.1   Introduction 

In  the  preceding  chapter,  the  CC  equations  were  derived  in  a  spin  orbital  framework. 
This  approach  allows  the  use  of  wide  range  of  reference  functions  (RHF,  UHF,  ROHF, 
QRHF,  etc.)  thus  allowing  application  to  both  closed  and  open  shell  systems  with  relative 
ease.  For  open  shell  molecules,  these  spin  orbital  based  CC  methods  may  be  characterized 
as  eigenfunctions  of  the  S2  operator  in  a  projected  sense,  (0  52  $M)  =  s(s  +  1),  for 
ROHF  and  QRHF  references,  or  not  at  all,  for  UHF  references.  Although  they  have  been 
very  successful  [34],  failure  to  be  a  rigorous  eigenfunction,  S2\$Cc)  =  s(s  -f  l)|\I>cc)> 
makes  this  approach  less  desirable  to  some  workers  in  the  field.  The  Fock  space  coupled 
cluster  method  grew  out  of  a  search  for  an  approach  to  open  shell  systems  which  would 
be  a  rigorous  eigenfunction  of  Ss  and  be  capable  of  handling  the  multiple  reference 
determinants  needed  to  represent  many  open  shell  states. 

There  are  two  primary  approaches  to  the  treatment  of  multireference  systems.  The 
Hilbert  space  approach,  also  referred  to  as  state-specific,  or  valance-universal,  treats 
a  manifold  of  states  with  a  constant  number  of  electrons,  and  each  Hilbert  space  CC 
calculation  focuses  on  a  particular  state.  The  Fock  space  (FS)  approach  treats  states  with 
different  numbers  of  electrons  (mathematically,  a  Fock  space  is  a  sum  of  Hilbert  spaces 
with  different  numbers  of  electrons),  simultaneously  calculating  several  different  states. 
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The  FSCC  method  is,  at  present,  the  more  developed  of  the  two,  and  offers  the  useful 
feature,  mentioned  in  Chapter  1,  that  it  can  provide  data  for  more  than  one  state  in  a 
single  calculation. 

In  this  chapter,  we  will  develop  the  FSCC  equations  in  more  detail,  obtaining  spin 
orbital  equations  for  IPs  at  the  FSCCSDT  level.  The  FSCC  approach  will  be  compared 
with  the  familiar,  and  (reasonably)  well  understood  single  reference  CC  approach  to 
point  up  the  inequivalence  between  FSCC  and  CC  calculations  at  the  same  level  of 
truncation  and  note  the  increased  importance  of  higher  excitations  (triples  in  particular) 
in  the  FSCC  approach.  Finally,  a  number  of  approximations  to  the  full  inclusion  of  triple 
excitations  are  proposed  in  order  to  include  some  triple  excitation  effects  but  reduce  the 
computational  effort  required. 

4.2  Structure  of  FSCC 

A  multireference  description  of  a  correlated  wavefunction  means  that  the  determi- 
nants expansion  is  partitioned  into  two  components.  The  model  space  includes  the 
determinants  desired  as  references,  and  the  complementary  space,  which  is  orthogonal 
to  the  model  space,  contains  excited  determinants  resulting  from  each  of  the  model 
space  functions.  The  same  division  holds  trivially  for  single  reference  CC,  with  the  one 
reference  determinant  being  the  model  space  and  all  excited  determinants  being  in  the 
complementary  space.  These  two  spaces  may  be  represented  by  projection  operators,  P, 
and  Q,  for  the  model  and  complementary  spaces  respectively,  defined  from  the  functions 
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they  include: 


V 


(4.1) 


which  are  orthogonal 


PQ  =  QP  =  0  (4.2) 


and  idempotent 


PP  =  P,  QQ  =  Q.  (4.3) 


For  convenience,  we  choose  a  single  determinant  one-electron  reference  function, 
|0),  to  provide  an  unambiguous  identification  of  occupied  and  virtual  orbitals.  The  idea 
of  a  particle-hole  rank,  denoted  (m,n),  is  used  to  indicate  the  states  related  to  |0)  by 
the  addition  of  m  electrons  [particles  in  the  otherwise  empty  virtual  orbitals)  and  the 
removal  of  n  electrons  (holes  in  the  occupied  orbitals).  Thus  singly  ionized  states  are 
in  the  (0,1)  sector,  electron-attached  states  are  in  the  (1,0)  sector,  and  the  (1,1)  sector 
includes  states  obtained  by  exciting  a  single  electron  from  |0).  Higher  particle-hole 
ranks  give  multiply  ionized  or  multiply  attached  states,  excitations  involving  more  than 
one  electron,  and  combinations  such  as  "shake-ups"  (combination  of  detachment  of  an 
electron  and  excitation  of  another).  The  idea  of  the  model  space  and  its  orthogonal 
complement  can  also  be  specialized  to  a  particular  particle-hole  rank,  giving  projection 
operators  p(m-n)  and  Q(m>n>,  respectively. 
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Figure  4.1.  Division  of  orbital  space  into  active  and  inactive  groups. 

In  order  to  define  the  model  space,  the  orbital  space  of  the  reference  determinant, 
|0),  is  divided  into  two  classes.  Active  orbitals  are  those  which  change  occupation  in  the 
various  model  functions,  while  inactive  orbitals  maintain  the  same  occupancy.  Since  the 
reference  function  has  already  been  divided  into  virtual  and  occupied  orbitals,  the  result 
is  a  four  different  categories  of  orbitals:  active  occupied,  active  virtual,  inactive  occupied, 
and  inactive  virtual.  These  divisions  are  shown  in  Figure  4.1.  Once  the  active  orbitals 
are  chosen,  the  model  space  for  each  particle-hole  rank  is  determined  —  all  appropriate 
occupations  of  the  active  space  are  included.  For  example  in  Figure  4.1,  the  (0,1)  sector 
(one  electron  removed)  model  space  includes  six  determinants  and  the  (1,0)  sector  has 
four  because  there  are  four  ways  to  add  a  single  electron  to  the  active  virtual  space; 
there  are  12  spin-conserving  single  excitations  within  the  active  space,  hence  12  model 
functions  in  the  (1,1)  sector.  It  should  be  noted  that  this  scheme  provides  for  exactly  one 
model  function,  |0),  in  the  (0,0)  sector,  and  the  calculation  in  this  sector  is  exactly  the 
single  reference  CC  method  discussed  in  the  preceding  chapter. 
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Model  spaces  can  be  divided  into  two  types,  complete  or  incomplete  [35].  A  complete 
model  space  is  one  that  contains  all  possible  occupations  of  the  active  orbitals,  while  in 
an  incomplete  model  space  (IMS),  some  occupations  of  the  active  orbitals  belong  to  the 
Q  space.  The  pure  ionization  and  pure  electron  addition  sectors,  (0,n)  and  (m,0)  sectors 
are  always  complete,  but  any  mixed  sector  is,  in  general,  incomplete.  For  example,  the 
(1,1)  sector  has  all  single  excitations  within  the  model  space  in  P,  but  places  all  double, 
triple,  and  higher  excitations  into  the  Q  space.  The  (1,1)  sector  is,  however,  part  of  a 
special  class  of  IMS  known  as  the  quasi-complete  model  space  [35,36].  In  this  case,  the 
active  orbitals  are  divided  into  disjoint  sets,  each  with  a  fixed  number  of  electrons.  If 
each  of  these  subsets  is  complete  within  itself,  the  entire  model  space  is  quasi-complete. 

The  FSCC  wavefunction  for  a  state  is  obtained  from  a  wave  operator,  ft,  operating 
on  a  linear  combination  of  model  space  functions: 

I^FSCC)  =  ni**>-  <4-4) 

where 

\*k)  =  52Cki\4>i).  (4.5) 

The  wave  operator  used  in  FSCC  is  an  exponential  of  cluster  operators,  just  as  in  the 
single  reference  theory,  but  with  a  normal  product  used  to  avoid  contractions  within  the 
cluster  operators  itself  (this  is  not  necessary  in  the  single  reference  formalism  since  it  is 
impossible  for  those  cluster  operators  to  contract  among  themselves), 


n 


{eS}.  (4.6) 
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The  cluster  operator  itself,  however,  has  a  somewhat  different  structure  from  the  single 
reference  case: 

§  —  g{m,n)  _     V^     y(*,0  (47) 

k=0.    m 
1  =  0.     n 

where  (m,n)  is  the  highest  sector  relevant  to  the  problem,  and  the  T^')  for  a  given  sector 
includes  single,  double,  and  higher  excitation  cluster  operators  for  that  sector: 

T(k,i)  =  T(k,l)  +  T(k,l)  +  T(k,i)  +  _  _  _  (4  g) 

Such  a  wave  operator  is  valance-universal  because  it  applies  to  all  sectors  from  (0,0)  to 
(m,n).  As  mentioned  above,  the  (0,0)  sector  is  a  CC  calculation  with  |0)  as  the  reference, 
so  r^0,0'  is  nothing  more  than  the  familiar  CC  cluster  operator  described  in  the  previous 
chapter.  Cluster  operators  for  other  sectors  are  different  from  the  (0,0)  ones  because  they 
apply  to  spaces  with  different  numbers  of  electrons,  or  at  least  electrons  in  different  places 
than  in  |0)  (i.e.  active  virtuals  may  contain  electrons  and  active  occupieds  may  receive 
electrons).  In  the  (0,1)  sector,  for  example,  one  of  the  virtual  orbital  labels  on  a  cluster 
amplitude  becomes  an  active  occupied,  corresponding  to  the  fact  that  in  this  sector,  an 
electron  has  been  removed  from  that  group,  making  it  possible  for  the  empty  orbital 
to  receive  an  electron.  The  differences  by  sector  in  the  cluster  operators  is  depicted  in 
graphical  form  by  Rittby  and  Bartlett  [37]  for  a  number  of  sectors,  and  some  implications 
of  these  differences  are  discussed  in  Section  4.3,  below. 

To  derive  the  FSCC  equations  themselves,  we  begin  from  the  Schrodinger  equation 
for  the  FSCC  wavefunction, 

Hn\yk)  =  Ek(2\yk).  (4.9) 
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This  equation  can  be  rewritten  in  terms  of  the  model  space  functions  by  defining  a 
matrix  t'  having  the  Et  as  eigenvalues  and  C  as  eigenvectors,  and  using  the  model  space 
projector,  P,  to  effectively  encompass  the  entire  model  space, 

HSlP  =  {lt'P.  (4.10) 

Using  the  fact  that  the  (0,0)  sector  cluster  operator  commutes  with  the  others  [38],  we 
can  decouple  the  (0,0)  sector  calculation  from  the  rest  of  the  FSCC  calculation.  First, 
we  extract  T'0'0'  from  the  wave  operator, 

„  r(0.0)-  -W0.0)  r     5(m,n)1 

ft  =  e1       ft  =  e1       ie*         >,  (4.11) 

where 

§(m,n)  _  J-  T{k'l],  (4.12) 

k=0...m,l  =  0      n 

and  multiply  both  sides  of  the  Schrodinger  equation  by  exp  f  —  T^0'0)]  to  obtain 

e-^      HNe^      OP  =  Ue'P  (4.13) 

Recalling  that  the  coupled  cluster  effective  Hamiltonian,  derived  in  the  preceding  chapter, 


is 

T(0.0)  t(0.°)  i      i    -   i      v 

HN  =  e~^      Her       —  <0|//  0> 


(4.14) 


T<0,0)  T<0,0) 


cc,\0) 


where  £Cc,|o)  1S  tne  CC  energy  of  the  reference  function  |0),  and  defining  a  new  matrix, 
e  =  e'  —  £CCi|o)'  which,  when  diagonalized,  represents  energy  differences  between  £CC||o) 
and  the  final  state,  we  can  write  the  FSCC  Schrodinger  equation  in  the  simplified  form: 

HN£lP  =  tteP.  (4.15) 
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Projecting  Eq.  4.15  on  the  left  with  P(l   l,  we  can  isolate  e  as  the  FSCC  effective 
Hamiltonian, 

PHNef£P  =  P&-1HNClP  =  PeP,  (4.16) 

which,  upon  diagonalization,  gives  the  state-to-state  energy  differences  we  seek.  Note  that 
this  result  could  have  been  obtained  by  requiring  intermediate  normalization,  PCtP  =  1, 
but  Mukherjee  [39]  has  shown  that  for  the  general  case  of  an  incomplete  model  space, 
intermediate  normalization  is  incompatible  with  size  extensivity.  Although  this  will  not 
cause  problems  for  the  (0,1)  sector,  our  ultimate  focus,  we  prefer  to  derive  generally 
applicable  equations  as  far  as  possible.  To  determine  the  cluster  operators,  we  project 
Eq.  4.15  on  the  left  with  the  complementary  space  determinants,  giving 


QHNttP-QttPHNetfP  =  0,  (4.17) 

which  is  often  referred  to  as  the  Bloch  equation.  Mukherjee  and  coworkers  have  shown 
that  for  complete  and  quasi-complete  model  spaces,  H N  eff  in  Eq.  4.17  is  connected  and 
the  energies  obtained  by  diagonalizing  it  are  size  extensive  [39-44]. 

Haque  and  Mukherjee  [38,45]  have  also  shown  that  there  is  a  subsystem  embedding 
condition  (SEC)  in  the  FSCC  equations  which  allows  Eq.  4.17  to  be  solved  separately 
for  each  sector  in  a  hierarchical  fashion.  As  already  observed,  the  (0,0)  sector  reference 
problem  decouples  from  the  rest  of  the  Fock-space  problem.  The  (1,0)  and  (0,1)  sectors 
can  be  solved  with  only  the  (0,0)  results  —  they  do  not  depend  on  any  higher  sectors. 
Given  the  (0,0),  (1,0),  and  (0,1)  sector  amplitudes  and  HNc^,  the  (2,0),  (1,1),  and  (0,2) 
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sectors  can  be  solved  —  once  again,  no  information  from  higher  sectors  is  required. 
In  this  fashion,  the  solution  of  the  desired  (m,n)  sector  is  built  up  from  lower  sectors, 
rather  than  requiring  the  solution  to  a  set  of  coupled  equations  for  all  sectors.  The  SEC 
is  also  very  useful  in  the  development  of  computer  programs  for  FSCC  because  the 
implementation  of  a  new  sector  should  not  require  any  changes  to  the  lower  sectors. 

4.3  Comparison  with  Single  Reference  CC  Theory 

Before  going  on  to  derive  the  more  specific  equations  for  Fock-space  methods,  it 
is  worth  taking  a  look  at  how  the  FSCC  approach  differs  from  the  more  familiar  single 
reference  formulation. 

As  has  been  noted,  the  Fock  space  coupled  cluster  amplitudes  differ  from  (0,0) 
sector  amplitudes  in  that  they  involve  excitations  between  the  occupied  or  virtual  orbitals 
and  orbitals  of  the  same  type  in  the  active  space.  In  the  reference  function,  |0),  such 
amplitudes  would  be  zero  because  active  occupied  orbitals  would  be  fully  occupied,  so 
no  excitation  could  deposit  an  electron  there,  and  active  virtuals  would  be  completely 
unoccupied,  providing  no  electron  to  excite  to  another  virtual.  In  the  model  space, 
however,  there  will  be  determinants  with  holes  in  the  active  occupied  space  and  particles 
in  the  active  virtual  space,  so  such  amplitudes  may  be  nonzero.  This  fact  changes  the 
character  of  FSCC  single,  double,  etc.,  excitation  amplitudes  from  their  familiar  (0,0) 
sector  interpretations  and  means  that  a  given  truncation  of  FSCC  will  not  provide  the 
same  results  for  a  given  state  as  the  same  truncation  of  the  (0,0)  sector  equations  applied 
to  the  model  space  function. 
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Figure  4.2.  Prototypical  reference  functions  and  model  space  determinants.     Active 
orbitals  are  shaded. 
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Figure  4.3.  Determinants  appearing  both  in  the  single  excitation  space  of  CCSD  calcu- 
lations on  the  four  model  functions  and  in  the  FSCCSD  singles  space. 

To  make  this  idea  more  concrete,  consider  a  FSCCSD  IP  calculation 
compared  to  CCSD  calculations  based  on  the  individual  determinants  in  the 
model  space.  Using  an  RHF  reference  function  with  three  occupied  orbitals 
and  one  virtual  orbital,  if  the  two  highest  occupied  orbitals  are  made  ac- 
tive,   we   have   four   determinants   in   the   model   space,    as   shown   in   Figure  4.2. 

A  CCSD  calculation  on  any  of  the  model  space  functions  has  seven  determinants  in 
the  single  excitation  space,  and  the  union  of  the  four  single  excitation  spaces  contains 
20  distinct  determinants.  FSCCSD,  on  the  other  hand,  restricts  single  excitations  to 
be  from  an  inactive  occupied  orbital  to  an  active  occupied,  which  allows  only  four 
amplitudes  and  yields  two  distinct  determinants,  shown  in  Figure  4.3.    The  remaining  18 


52 


T 
1 

T 
T  i 

T 
T 

T 

T  i 

T  I 

T  i 

T  1 

I 

T  1 

T 

t24 

l23 

,34 
l32 

,24 
r21 

t24 

22 

t24 

Z34 

l32 

r31 

Figure  4.4.  Examples  of  determinants  arising  from  spectator  double  excitations  in  FSCC, 
corresponding  to  single  excitations  in  the  CC  model.     Shown  are  all  determinants 
generated  from  the  |^>j  '    )  model  function  with  the  relevant  amplitudes  from  T2l  ' 
indicated  below  each  one.  Orbitals  are  numbered  from  lowest  energy  (bottom)  to  highest 
(top);  beta  spin  electrons  are  indicated  by  a  bar  over  the  orbital  index. 

CCSD  determinants  are  generated  in  the  Fock  space  calculation  by  a  class  of  excitations 
referred  to  as  spectator  doubles.  FS  double  amplitudes  include  one  excitation  from  an 
occupied  orbital  to  an  active  occupied.  There  is  no  restriction  on  the  initial  orbital,  so 
it  may  be  active  as  well  —  in  fact  it  may  be  the  same  as  the  final  orbital.  This  null 
excitation,  which  is  necessary  to  distinguish  between  excitations  connecting  different 
model  functions  to  the  same  final  determinant  [46,45],  makes  spectator  doubles  effectively 
single  excitations  in  the  usual  coupled  cluster  sense.  Figure  4.4  shows  the  spectator 
doubles  arising  from  model  function  \(f>\  '    ),  along  with  the  corresponding  amplitudes. 

When  considered  at  truncation  to  only  double  excitations  (CCD  and  FSCCD),  the 
Fock  space  approach  is  going  to  include  some  (but  not  all)  effects  of  single  excitations  — 
a  potential  advantage.  When  higher  cluster  operators  are  included,  however,  FSCC  is 
at  an  unequivocal  disadvantage.  Compare,  for  the  same  model  problem,  the  double  ex- 
citations obtained  from  the  CC  and  FSCC  approaches.    FSCCSD  generates  22  "true" 
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Figure  4.5.  Examples  of  determinants  arising  from  spectator  triple  excitations  in  FSCC, 
corresponding  to  double  excitations  in  the  CC  model.     Shown  are  all  determinants 
generated  from  the  \(f>\  '    )  model  function  with  the  relevant  amplitudes  from  T3  ' 
indicated  below  each  one.  Orbitals  are  numbered  from  lowest  energy  (bottom)  to  highest 
(top);  beta  spin  electrons  are  indicated  by  a  bar  over  the  orbital  index. 

doubly  excited  determinants,  while  CCSD  calculations  on  the  individual  model  func- 
tions yield  18  additional  determinants.  The  "extra"  CCSD  determinants  are  created 
in  FSCC  by  spectator  triple  excitations,  a  few  of  which  are  displayed  in  Figure  4.5. 


From  these  examples,  it  should  be  evident  that  for  a  FSCC  IP  calculation  to  match 
the  expansion  space  of  a  single  reference  CC  of  a  given  truncation  on  one  of  the  model 
functions,  the  FSCC  calculation  must  include  spectator  contributions  from  the  next  higher 
cluster  operator.  The  same  would  be  true  of  EA  calculations,  and  in  higher  sectors,  more 
indices  of  a  given  amplitude  terminate  in  the  active  space  (p=m+n  for  the  (m,n)  sector), 
so  it  would  be  necessary  to  include  spectators  from  the  p^  higher  cluster  operator  to 
obtain  an  expansion  space  equivalent  to  that  of  the  single  reference  CC. 
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4.4  Basic  Equations 

From  this  point,  we  will  focus  on  the  (0,1)  sector,  which  describes  one-electron 
ionizations  with  respect  to  the  reference  state  and  are  of  particular  interest  in  this  work. 

From  Eq.  4.17,  the  detailed  operator  equations  which  determine  each  cluster  operator 
are  easily  determined  by  expanding  the  exponential  wave  operator,  as  in  the  previous 
chapter.  Note  that  in  FSCC,  products  of  cluster  operators  in  the  exponential  expansion 
contribute  to  other  sectors,  making  the  (0,1)  equations  a  good  deal  simpler  than  their 
CC  counterparts. 


^0'1)^,eff^(°'1)  =  i,(0'1){/jv(l  +  r{0,1)  +  ^0,1)) 

+^(^f1)  +  rf1))+/7^rj0'1)}cJp(°'1) 

o  =  Q(W{fN(i  +  rr]  +  if1')  +  wn(t^  +  if1') 
+IIINT^  -  T^R^}pm 


(4.18) 


(4.19) 


o  =  Q<tu)  {fN  (t2(0,1)  +  if 1})  +  wN  (i  +  if  ^  +  rt1]  +  ?f 1}) 


(4.20) 


55 

+  /7/iV(i  +  r1(0'1)  +  ^0'1)  +  r3(0'1))  (4.2D 

_r(o,i)i/(o,i)  1    p(0,l) 
3         iv.eff  J  c 

4.5  Reduction  to  Spin  Orbital  Equations 

Equations  4. 18—4.21  are  transformed  to  spin  orbital  form  following  the  same  proce- 
dures used  in  the  previous  chapter  for  the  CC  amplitude  equations  and  H^.  The  resulting 
equations  for  the  FSCCSDT  method  for  ionization  potentials  are  given  in  Eq.  4.22—4.25, 
below.  The  notation  is  the  same  as  that  used  in  the  previous  chapter  with  a  few  additions. 
The  Fock-space  cluster  amplitudes  are  called  r  to  avoid  any  possibility  of  confusion  with 
the  (0,0)  sector  amplitudes,  t.  A  single  dot  on  an  orbital  index,  m,  indicates  restriction  to 
inactive  orbitals,  and  a  double  dot,  m,  indicates  restriction  to  active  orbitals  (this  choice 
of  notation  is  meant  to  remind  readers  familiar  with  the  diagrammatic  notation  of  the 
"circled  arrow"  used  to  denote  inactive  lines,  and  the  double  arrow  for  active  lines). 


-'  f      Je  1  /„,  ,il;"   \Je  *■/  f\Jef 


HN,e«)ijc  =  fa  ~  4iT*  +  /"•*£.  -  oi^^Kn  +  l{mn\\ef)T^n         (4.22) 


0  =  4-4^I  +  ^rt-i(m,^e)ri  +  i(mnHe/)r^n  +  rf  {HN^   (4.23) 
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0  =  Ua\\ki)  +  faer^  -  (-lYP(ki)fmiriam  -  (ma||M)rj 


+  ^(-iyP(ki)(mn\\ki)rtn  -  (ma\\ke)rti  -  l-{man\\kze)r^n 

+  fmertn  +  \{-lf  P{ki)(mn\{ke)T}Zi  +  \(ma\\ef)Tti  +  rg"  (^.df)?™ 

(4.24) 


0  =  (lab\\krj)  -  (-)pP(k/ij)(mb\\ij)Tlfm  +  (-)pP(*z/j)(a%J>£ 
+  l-(-)pP(ab)(m.an\\kij)T±  -  (-)pP(ki/ < j){mib\\kie)T%j  -  (mab\\kij)ri 
+  {-)pP(ab)fheTiTi  ~  i-)PP{kilj)fmAaL  +  (-)pP(ki/j)(mn\\ki)rtj 


(-)pP(^';,a6)<ma|>e>T^  +  T™b{HN^ 


lm 

(4.25) 


4.6  Approximations  to  the  Full  FSCCSDT  Model 

In  the  single  reference  CC  world,  the  full  CCSDT  model  has  seen  limited  use  for 
a  number  of  reasons 

•  Complexity  of  the  method  —  the  triples  equations  have  many  more  terms  than  single 
and  double  excitations,  thus  requiring  substantially  more  effort  to  code  them  into  a 
computer  program. 

•  Primary  and  secondary  memory  (i.e.  physical  memory  and  disk  space)  requirements 
for  storage  and  manipulation  of  the  cluster  amplitudes,  which  are  significantly  more 
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numerous  than  the  double  excitation  amplitudes  —  a  factor  of  n(N-n)  in  the  language 
of  section  3.6. 

Time  constraints,  since  many  more  amplitudes  must  be  determined,  the  time  required 
increases  in  a  similar  fashion. 

As  a  result,  a  number  of  approximations  to  the  full  CCSDT  method  have  been  proposed 
and  have  proven  quite  effective  at  incorporating  important  effects  of  triple  excitations 
while  at  the  same  time  minimizing  the  above  concerns  [33,47-52].  These  concerns 
might  also  be  raised  with  respect  to  the  full  FSCCSDT  method,  so  in  the  same  spirit, 
we  propose  several  approximations  of  FSCCSDT  as  a  compromise  between  the  cost  and 
performance  of  the  full  model. 

There  are  two  different  approaches  to  approximating  FSCC  triples  which  may  be 
considered.  One,  in  analogy  to  the  approximate  triples  treatments  of  single  reference  CC, 
involves  truncations  and  perturbative  approximations  to  the  triples  equations.  The  second 
method  is  to  restrict  the  triples  amplitudes  which  are  evaluated  to  something  less  than 
the  full  set.  This  approach,  which  has  not  been  employed  in  CC  theory,  is  suggested  by 
the  earlier  observation  that  spectator  triples  must  be  included  to,  in  some  sense,  duplicate 
the  expansion  space  of  a  single  reference  CCSD  calculation.  The  two  approaches  are  not 
exactly  orthogonal,  but  can  be  considered,  by  and  large,  independently. 

In  single  reference  CC  theory,  there  are  no  self-evident  subsets  of  the  triples 
amplitudes,  so  this  approach  has  not  been  tried  (or  at  least  not  reported  in  the  literature). 
In  FSCC,  however,  there  are  two  obvious  subsets  which  might  be  considered.  The 
smallest  set  would  be  the  spectator  amplitudes,  as  suggested  by  the  discussion  in  the 
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previous  section.  A  larger  set  would  be  to  include  all  active-active  amplitudes.  This  set, 
which  includes  the  spectators,  allows  excitations  from  any  active  orbital  into  any  other 
(spectators  involve  an  active  orbital  scattering  into  itself,  and  the  full  set  of  triples  allows 
scattering  from  inactive  orbitals  to  occupied  orbitals  as  well  as  active  to  active). 

As  we  saw  in  the  previous  section,  it  is  necessary  to  introduce  spectator  triples  in  the 
Fock  space  approach  in  order  to  match  the  expansion  space  of  a  CCSD  treatment  of  each 
of  the  model  functions  in  turn.  On  the  other  hand,  since  FSCC  describes  all  states  of 
interest  with  a  single  set  of  amplitudes  (with  diagonalization  of  the  effective  Hamiltonian 
distinguishing  among  individual  states),  it  can  be  argued  that  it  is  impossible  for  it  to 
be  equivalent  to  a  separate  CC  calculation  for  each  model  function,  and  the  differences 
in  the  expansion  space  between  the  two  approaches  should  not  be  a  concern.  Indeed,  a 
number  of  calculations  in  the  literature  [53,54,43,55,  for  example]  support  the  idea  that 
FSCCSD  is  well-balanced  and  the  results  obtained  are  numerically  on  a  par  with  other 
methods.  The  CCSD  method  has  been  widely  applied,  and  its  behavior  is  reasonably  well 
understood  —  certainly  there  is  far  more  experience  with  CCSD  than  with  any  FSCC 
method  —  and  this,  more  than  any  special  quality  of  CCSD  underlies  the  desire  for  a 
FSCC  method  that  can  mimic  it.  Only  a  critical  examination  of  the  proposed  spectator 
and  active-active  approximations  for  a  wide  range  of  systems  can  determine  if  they  will 
improve  the  results  over  FSCCSD  or  destroy  the  balance  of  the  FSCC  method. 

The  second  approach  to  reducing  the  cost  of  a  FSCCSDT  calculation  is  to  approximate 
or  ignore  the  most  expensive  terms  to  evaluate.  In  this  sense,  the  simplest  possible 
approximation  would  be  to  evaluate  the  triple  excitation  contribution  perturbatively  from 
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converged  FSCCSD  amplitudes.    The  single  reference  analog  to  this  approach  is  the 
CCSD+T(CCSD)  method,  in  which  the  triple  excitation  contribution  is  correct  through 
fourth  order  —  the  first  order  at  which  it  appears. 

For  the  FSCC  perturbation  expansion,  we  follow  the  choice  of  Pal,  Rittby,  and  Bartlett 
[56]  which  was  also  implicit  in  the  work  of  Haque  and  Kaldor  [57],  of  counting  orders 
in  terms  of  H,  the  reference  state  Hamiltonian.  Alternatives,  such  as  using  H  would  lead 
to  inconsistencies  between  the  (0,0)  sector  definitions  and  those  of  higher  sectors.  For 
example,  /  and  W  are  first  order  in  //,  while  the  higher  components,  777, ...  are  at  least 
second  order.  Counting  orders  in  H  would  place  all  of  these  terms  on  the  same  level, 
making  a  perturbation  expansion  of  Eq.   4.21  useless. 

Although  in  both  CC  and  FSCC  the  triple  excitation  amplitudes  first  appear  in  second 
order  of  perturbation  theory,  in  the  CC  formalism  excitations  beyond  singles  and  doubles 
contribute  only  indirectly  to  the  energy  via  their  contributions  to  the  singles  and  doubles 
amplitudes,  and  thus  triples  first  appear  in  the  energy  in  fourth  order.  In  FSCC,  all  clusters 
contribute  directly  to  HN  ejy  and  thus  the  energy,  so  the  triples  first  manifest  themselves 
in  third  order.  Inclusion  of  the  third  order  triples  effects  based  on  converged  FSCCSD 
amplitudes  was  first  proposed  by  Haque  and  Kaldor  [57],  and  called  CCSD+T.  Later,  Pal, 
Rittby,  and  Bartlett  [56]  discuss  the  third  order  model  using  the  name  MRCCSD+T(3) 
and  propose  several  other  approximations.  The  MRCCSD+T*(3)  model  attempts  to 
include  some  higher  order  terms  by  replacing  the  /  and  W  required  in  T(3)  with  their 
counterparts  from  the  effective  Hamiltonian,  /  and  W.  The  nature  of  this  method,  which 
will  be  referred  to  as  FSCCSD+T(3)  in  this  work,  can  be  seen  from  the  spin  orbital 


60 
FSCCSDT  equations  given  above.  The  terms  which  contribute  to  the  TJ  '  '  amplitudes 
in  second  order  are  those  on  the  first  row  of  Eq.  4.25,  with  appropriate  restrictions  on 
the  Htf  elements:  the  W^  elements  of  the  second  and  third  terms  are  simply  the  bare 
two-electron  integrals  in  first  order,  and  only  two  terms  of  (lab\\kij)  appear  at  this  level 
(the  first  two  terms  of  (iab\\jkl\  in  Table  3.3).  The  second  order  triples  amplitudes  in 
turn  contribute  to  //^,eff  (and  thus  the  energy)  in  third  order  thought  the  final  term  of 
tf^eff  in  Eq.    4.22. 

The  extension  of  this  approach  to  higher  orders  of  perturbation  theory  is  obvious.  A 
fourth-order  method,  FSCCSD+T(4),  would  be  obtained  by  the  following  steps: 

1.  Second-order  T3  '  amplitudes,  as  evaluated  above,  give  third-order  contributions 
when  inserted  in  the  Tj  '  and  T^  '  amplitude  equations.  This  correction  to  the 
Tj  '  and  r2  '  amplitudes  can  be  directly  contracted  into  H^^,  giving  a  fourth 
order  correction  to  the  energy. 

2.  The  third-order  T\  '  must  be  evaluated.  This  requires  the  fourth,  fifth  and  sixth 
terms  from  Eq.  4.25,  and  additional  terms  from  (/a6||A:z;)  (terms  3-10  of  (iab\\jkl) 
in  Table  3.3).  Note  that  even  if  the  (0,0)  reference  calculation  does  not  include  triple 
excitations,  approximate  T$  '  amplitudes  can  be  formed  at  second  order  in  pertur- 
bation theory,  and  therefore  must  be  included  in  the  general  case.  Computational 
experiments  will  be  required  to  determine  the  effect  of  approximate  T3  '     in  H^. 

Of  course  such  a  perturbation  expansion  could  go  on  for  a  very  long  time,  but 
in  practice,  better  results  are  likely  to  be  obtained  (based  on  SRCC  experience)  by 
iterating  an  approximate  triples  equation  as  part  of  the  FSCCSD  calculation.     The 
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simplest  approximation  for  which  this  makes  any  sense  is  T(4),  described  in  the  preceding 

paragraph  —  it  is  the  first  perturbative  method  for  which  T3       — >  T2  '     and  T3       — * 

T{  '     contributions  occur.    The  procedure  is  simple:    In  each  iteration  of  the  FSCC 

calculation,  after  the  single  and  double  excitation  amplitudes  have  been  evaluated,  the 

T(4)  procedure  is  used  to  calculate  the  contribution  of  triples  excitations  to  H^^  and 

the  singles  and  doubles  amplitudes.    Instead  of  this  being  the  end  of  the  calculation, 

the  new  TJ  '  '  and  T2(  '     amplitudes,  and  //jv,eff  are  fed  back  into  the  next  cycle  of 

iteration.    Note  that  there  are  no  T3        — >  T3  '     contributions  in  this  method,  so  the 

triple  excitation  amplitudes  do  not  need  to  be  kept  from  one  iteration  to  the  next,  and 

also  avoids  the  most  expensive  part  of  the  T3  '    equation  to  evaluate.  This  approach  is  the 

same  as  that  taken  by  Urban  et  al.  [48]  in  approximating  single  reference  CCSDT.  In  the 

hopes  of  maintaining  some  consistency  in  nomenclature  between  SRCC  and  Fock-space 

approaches,  this  method  should  be  called  FSCCSDT-1. 

In  the  single  reference  case,  the  idea  has  been  extended  in  a  series  of  additional 
iterative  approximations,  CCSDT-2  [48,50],  CCSDT-3  [48,50],  and  CCSDT-4  [51], 
which  include  more  and  more  of  the  terms  due  to  triples,  but  continuing  to  exclude 
r3  — >  T3  ,  thus  avoiding  the  most  expensive  terms  both  computationally  and  for 
storage  space.  In  the  FSCC  case,  the  T3  '  '  equation  has  a  simpler  structure,  and  only 
Tj  '     — *  T3  '     terms  remain  after  the  FSCCSDT-1  approximation. 

Over  the  years,  experience  with  various  SRCC  approximate  triples  methods  has 
led  to  refinements  in  several  of  them.  One  of  the  main  observations  is  that  in- 
cluding selected  higher  order  terms  in  a  perturbative  or  iterative  method  can  pro- 
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duce  better  results.  For  example,  Watts,  Gauss,  and  Bartlett  [58]  have  recently  intro- 
duced the  CCSD[T]  method,  which  is  a  modification  of  CCSD+T(CCSD),  the  sim- 
plest SRCC  perturbative  triples  method.  They  differ  in  the  inclusion  of  a  fifth- 
order  term  and  in  the  requirement  that  non-Hartree-Fock  (i.e.  QRHF  or  ROHF)  ref- 
erences be  treated  properly.  This  requires  that  terms  involving  off-diagonal  /  ele- 
ments be  considered  since  they  are  zero  for  a  Hartree-Fock  reference,  but  nonzero 
in  the  non-HF  case.  Transforming  the  reference  functions  to  semicanonical  orbitals 
[59]  has  the  effect  of  diagonalizing  the  occupied-occupied  and  virtual-virtual  blocks 
of  the  Fock  matrix,  so  that  terms  involving  ftJ  or  faf,  can  be  avoided.  This  leaves 
occupied-virtual  Fock  matrix  elements,  fat  to  be  accounted  for,  which  involves  shift- 
ing these  terms  down  one  order  of  perturbation  theory,  from  first  to  zeroth  or- 
der. 

In  the  Fock-space  coupled  cluster  framework,  however,  the  situation  is  somewhat 
different.  The  transformation  to  semicanonical  orbitals  diagonalizes  blocks  of/,  not  fa. 
Even  with  such  a  transformation,  the  leading  terms  in  fa  and  fab  are  first  order  in  the  non- 
HF  case.  This  would  result  in  FT^  '  — >  T3  contributions,  one  of  the  most  expensive 
terms  in  FSCCSDT,  beginning  with  the  third-order  amplitudes  (FSCCSD+T(4)).  Thus 
as  far  as  approximations  to  FSCCSDT  go,  such  a  requirement  would  be  reasonable  only 
for  the  most  basic  third-order  method,  FSCCSD+T(3).  This  would  add  the  FT^'Q)  term 
of  (lab\\kij)  (see  Table  3.3)  and  the  FT2(0'0)  terms  of  (mb\\ij)  and  (ab\\ej)  (see  Table 
3.2),  and  would  require  evaluation  of  the  FT^0,1)  ->  r2(0,1)  -»  H(°'^  term  of  Eq.  4.24 
since  fme  is  now  zeroth  order. 
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There  is  one  further  matter  regarding  approximation  to  the  FSCCSDT  equations. 
Mattie  [60]  has  proven  that  for  the  (0,1)  sector  of  Fock-space,  the  values  of  the  roots 
obtained  by  diagonalizing  ///v,eff  are  independent  of  the  active  space  used  for  the 
calculation.  In  other  words,  if  two  orbitals,  a  and  b  are  taken  as  active,  the  resulting 
ionization  potentials  would  be  identical  to  those  obtained  from  separate  calculations  with 
either  a  or  b  alone  active.  This  is  a  very  useful  result  because  it  means  that  we  do  not 
have  to  worry  about  choosing  the  "right"  active  space  for  a  given  calculation  in  order 
to  get  good  results.  However  the  proof  of  this  invariance  rests  on  the  FSCC  amplitudes 
and  i/^,eff  satisfying  the  Fock-space  Bloch  equation  4.17.  None  of  the  approximate 
FSCCSDT  methods  described  above  satisfy  Eq.  4.17,  so  the  invariance  is  lost.  It  remains 
to  be  seen  how  large  this  effect  will  be,  or  if  there  is  an  alternative  route  to  the  invariance 
condition  that  will  work  with  approximate  solutions. 

As  a  start  towards  understanding  the  FSCC  triple  excitation  correction  as  well  as  we 
now  understand  the  single  reference  triples  correction,  the  remainder  of  this  work  will 
focus  on  the  third-order  triples  approximations,  FSCCSD+T(3)  and  FSCCSD+T*(3).  The 
implementation  of  these  methods  for  both  closed  and  open-shell  reference  functions  is 
described  in  some  detail  in  the  Appendix,  including  the  final,  spin-integrated  equations 
from  which  the  computer  code  was  written.  In  the  following  chapters,  these  methods  are 
applied  to  a  variety  of  molecules  in  order  to  assess  their  performance. 


CHAPTER  5 

A  CRITICAL  ASSESSMENT  OF  FSCCSD  AND 

THE  THIRD-ORDER  TRIPLES  APPROXIMATIONS 


5.1   Introduction 

Before  applying  a  new  theoretical  method  to  otherwise  unknown  molecular  systems,  it 
is  important  to  get  a  feel  for  its  behavior  —  strengths,  weaknesses,  and  overall  quality  — 
by  comparing  the  results  obtained  with  new  method  to  other  theoretical  and  experimental 
data  for  similar  types  of  problems.  This  "calibration"  of  a  method  is  a  continuous  and 
dynamic  process  —  each  time  the  method  is  used  for  a  new  problem,  more  information 
about  its  performance  is  gained,  thus  modifying  the  assessment  of  its  accuracy  and 
applicability. 

Choosing  the  appropriate  standard  of  comparison  for  a  new  method  is  in  most  cases 
not  a  simple  matter.  Experimental  results  are  an  obvious  choice,  but  there  some  potential 
problems.  Experimental  corresponding  to  readily  computed  results  are  not  always 
available  or  may  not  be  known  accurately.  For  example,  FSCC  most  readily  provides 
vertical  ionization  potentials,  but  only  adiabatic  may  be  available  from  experiment.  For 
potential  energy  surfaces  —  geometries,  relative  energies,  and  the  like  —  there  is  usually 
very  limited  experimental  data.  Structures  may  be  known  for  minima  on  the  PES,  but 
experiments  which  can  probe  transition  states  are  quite  new,  and  have  not  yet  produced 
a  sizeable  body  of  data  on  different  systems.  The  vibrational  spectra  of  a  given  species 

64 


65 

may  be  only  partially  known,  and  in  many  experiments,  assignment  of  spectral  lines  to 

particular  species  can  be  a  very  hard  problem.  It  is  also  possible  for  a  calculation  to  obtain 
the  right  (experimental)  answer  but  for  the  wrong  reasons  —  a  fortuitous  cancellation  of 
errors.  If  such  problems  are  given  due  attention,  experimental  results  can  be  a  very 
useful  measure  the  performance  of  a  theoretical  method,  but  it  is  worth  considering  other 
standards  as  well. 

Comparison  of  a  new  method  with  experimental  results  often  introduces  more  vari- 
ables than  one  would  like.  For  example,  the  choice  of  basis  set  may  not  be  adequate  to 
obtain  high-quality  theoretical  results.  Full  CI  provides  the  exact  answer  for  a  given  basis 
set  and  Hamiltonian,  and  therefore  offers  an  excellent  standard  of  comparison  for  new 
methods.  Unfortunately,  FCI  calculations  are  very  expensive,  so  the  results  available  for 
comparison  are  somewhat  limited.  Another  option  is  to  compare  with  methods  that,  from 
a  well-defined  theoretical  standpoint,  are  better  than  the  method  in  question.  An  example 
of  this  would  be  comparing  methods  which  approximates  triple  excitation  effects  to  the 
full  triples  model  such  as  CCSD(T)  vs.  CCSDT,  or  FSCCSD+T(3)  vs.  FSCCSDT.  In 
this  work,  however,  as  is  often  the  case,  the  approximate  method  has  been  implemented 
before  the  full  method  because  of  the  additional  complexity  of  the  full  model,  forcing 
such  comparisons  to  wait  until  the  more  complete  method  is  implemented.  A  final  al- 
ternative is  comparison  of  the  new  method  with  one  for  which  the  behavior  is  relatively 
well  understood,  but  which  cannot  be  clearly  designated  as  a  "better"  method.  Although 
these  comparisons  are  the  least  definite,  they  are,  of  necessity,  the  most  common  in  the 
calibration  of  any  new  method. 


66 

This  chapter  is  a  beginning  at  the  calibration  of  the  FSCCSD+T(3)  and 
FSCCSD+T*(3),  models.  The  molecules  chosen  for  these  calculations  are  drawn  largely 
from  applications  of  the  FSCCSD  method  already  in  the  literature,  and  in  this  work 
extended  to  include  triple  excitation  contributions.  The  intent  is  to  get  an  overview  of 
the  performance  of  the  new  methods  for  a  variety  of  systems,  most  of  which  have  also 
been  the  subject  of  single  reference  CC  calculations,  thus  providing  a  useful  point  of 
comparison. 

5.2  Effect  of  Changes  in  the  Active  Space 

The  triples  approximations,  as  was  noted  in  the  previous  chapter,  do  not  have 
the  same  invariance  with  respect  to  the  choice  of  active  space  as  occurs  for  com- 
plete methods.  In  order  to  simplify  later  comparisons,  we  first  examine  the  effect  of 
the  change  in  active  space  on  the  third-order  triples  models.  The  "error  bars"  ob- 
tained in  this  section  may  then  be  kept  in  mind  later  on,  where  a  single  choice  of 
active  space  will  be  used  in  comparisons  of  the  new  FSCC  methods  with  other  re- 
sults. 

Because  //jv)eff  is  block  diagonal  due  to  molecular  symmetry,  the  active  space  within 
a  particular  symmetry  block  can  be  changed  without  affecting  the  others.  If,  in  a  series 
of  calculations,  the  active  space  for  a  particular  irrep  remains  unchanged,  the  IPs  of  that 
symmetry  will  not  change.  This  is  the  case  regardless  of  changes  to  the  active  space 
in  other  irreps.  Conversely,  changing  the  active  space  within  a  given  irrep  will  result 
in  slightly  different  IPs. 
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In  order  to  assess  the  magnitude  of  this  effect,  calculations  have  been  carried  out  on 
a  number  of  molecules  using  different  active  spaces.  The  molecules,  listed  in  Table  5.1, 
include  a  number  for  which  the  IPs  themselves  are  studied  below,  and  several  additional 
systems  with  fairly  low  symmetry  actual  or  computational  symmetry.  The  idea  behind 
these  additional  molecules  was  to  see  if  having  more  occupied  orbitals  in  each  irrep  (and 
therefore  a  larger  set  of  possible  active  spaces)  would  result  in  larger  variations  in  the 
calculated  IPs.  FSCCSD+triples  calculations  were  carried  out  on  each  molecule  using  a 
DZP  basis  at  the  DZP/SCF  optimized  geometry  for  a  set  of  active  spaces  consisting  of 

1.  the  largest  active  space  for  which  the  underlying  FSCCSD  could  be  converged 
(referred  to  as  the  "reference"  active  space), 

2.  single  orbital  active  spaces  for  each  of  the  orbitals  in  the  reference  active  space,  and 

3.  a  series  of  active  spaces  obtained  by  starting  with  the  lowest  IP  and  successively 
adding  the  next  highest  IP. 

This  choice  provides  a  slightly  larger  set  of  data  on  the  variation  of  lower  IPs,  but  those  are 
usually  the  roots  of  greatest  interest  anyway.  Analysis  of  these  calculations  involved  first 
reducing  the  data  to  obtain  a  set  of  unique  active  spaces  within  each  symmetry  block, 
then  computing  for  each  IP,  the  mean,  range,  and  mean  difference  from  the  available 
data.  The  mean  was  chosen  as  the  reference  point  because  there  is  no  theoretical  basis 
for  preferring  one  active  space  over  another  as  the  "correct"  one,  except  perhaps  for 
having  all  occupied  orbitals  active,  which  proved  impossible  to  converge  in  most  cases. 
Table  5.1  presents  the  largest  range  and  mean  difference  observed  among  the  various  IPs 
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Table  5.1.  Effect  of  changes  in  active  space  in  the  T(3)  and  T*(3)  models,  for  a  variety 
of  molecules.  Ranges  and  mean  differences  of  IPs  are  given  in  mH. 


Molecule 

Point 
groupa) 

Active 
orbitalsb) 

nc> 

FSCCSD+T(3) 

Mean 
Ranged) 

Diff.e) 

FSCCSD+T*(3) 

Mean 

Ranged) 

Diff.e) 

CaH^ 

cs 

3 

3 

0.023 

0.010 

0.023 

0.010 

CH2g) 

C2v 

3 

3 

0.078 

0.035 

0.065 

0.029 

CH2NH^ 

Cs 

4 

4 

0.744 

0.220 

0.545 

0.165 

CH2PH8' 

Cs 

4 

2/4 

0.103 

0.052 

0.074 

0.028 

CH3Clf) 

Cs  (C3v) 

4 

4 

0.008 

0.004 

0.006 

0.003 

CH3CNn 

Cs  (C3v) 

4 

4 

0.354 

0.133 

0.284 

0.107 

CH3NCf) 

Cs  (C3v) 

4 

4 

0.044 

0.022 

0.044 

0.022 

CH3NH20 

Cs 

4 

4 

0.226 

0.083 

0.175 

0.066 

H2C08> 

C2v 

2 

2 

1.118 

0.559 

0.931 

0.466 

HOC10 

Cs 

3 

2 

1.844 

0.921 

1.581 

0.790 

HOFn 

Cs 

3 

2 

3.467 

1.734 

2.757 

1.379 

N2e) 

D2h  (Pooh) 

3 

3 

1.297 

0.576 

0.861 

0.383 

a)  Computational  point  group  (full  point  group). 

b)  Largest  active  space  used  in  the  irrep  containing  the  root  showing  the  largest  variation. 

c)  Number  of  different  active  spaces  sampled  for  the  root  showing  the  largest  variation. 
Same  for  T(3)/T*(3)  unless  two  values  are  given. 

d)  Largest  range  (difference  between  largest  and  smallest  values)  for  a  single  root,   (in 
mH). 

e>  Mean  difference  (~  Y^=i  |^>  —  IF*  | )  for  the  root  showing  the  largest  variation  (in  mH). 
"  Calculated  with  a  DZP  basis  at  the  DZP/SCF  optimized  geometry. 
g)  Calculated  using  the  geometry  and  basis  used  in  the  ionization  potential  calculations 
in  the  next  section. 


for  each  molecule,  along  with  an  indication  of  the  number  of  distinct  active  spaces  in 
which  that  IP  was  computed. 
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The  largest  variation  is  observed  for  HOF,  with  HOC1,  N2,  and  H2CO  close  behind. 
For  a  number  of  other  species,  however,  the  variation  is  quite  small.  Interestingly,  the 
largest  variations  are  seen  among  molecules  with  higher  symmetry  and  smaller  active 
spaces  rather  than  large  ones.  This  obviously  suggests  that  the  number  of  occupied  or- 
bitals  is  not  the  most  important  determinant  of  the  variation,  however  examination  of  the 
FSCCSD  and  reference  CCSD  results  provided  no  clues  to  distinguish  molecules  like 
HOF  or  H2CO  from  ones  with  a  small  variation,  such  as  CH3CI. 

Comparison  of  the  T(3)  and  T*(3)  results  shows  that  the  range  and  mean  difference 
for  T*(3),  which  includes  some  higher-order  terms,  are  at  least  as  good  as,  and  generally 
better  than,  those  for  T(3).  Although  there  are  a  few  cases  where  going  to  the  T*(3) 
approximation  changes  which  of  the  IPs  has  the  largest  variation,  generally  it  is  the  same 
ionization  for  both. 

In  order  to  verify  that  the  variation  is  not  due  to  an  inadequacy  of  the  DZP  basis  sets 
used  for  most  results  in  Table  5.1,  a  second  set  of  calculations  was  carried  out  on  HOF 
using  Dunning's  correlation  consistent  PVQZ  basis  [61],  which  is  5s4p3d2flg  on  O  and  F, 
and  4s3p2dlf  on  H  for  a  total  of  125  basis  functions  (using  spherical  harmonic  functions), 
compared  to  4s2pld  and  2slp  for  the  DZP  basis,  totalling  37  basis  functions  (using 
cartesian  functions).  Table  5.2  presents  the  mean  ionization  potential  and  its  variation 
with  active  space  for  the  lowest  five  IPs  of  HOF  calculated  using  both  basis  sets.  There  are 
changes  in  the  range  of  variation  from  one  basis  to  the  other  —  some  roots  have  a  larger 
range  with  the  PVQZ  basis,  others  have  a  smaller  range.  In  general,  however,  the  they  are 
roughly  the  same,  suggesting  that  the  variation  is  largely  independent  of  the  basis  set. 
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Table  5.2.  Effect  of  changes  in  the  active  space  on  the  calculated  ionization  potentials 
(in  mH)  of  HOF  for  two  basis  sets. 


DZP  basis 

PVQZ  basis 

Sym. 

n 

Mean 
IP 

Range 

Mean 
Diff. 

Mean 
IP 

Range 

Mean 
Diff. 

FSCCSD+T(3) 

a' 

3 

545.7 

1.158 

0.608 

559.0 

1.176 

0.480 

a' 

3 

672.8 

3.186 

1.064 

686.0 

3.529 

1.281 

a' 

2 

775.7 

3.467 

1.734 

765.3 

3.714 

1.856 

a" 

2 

466.5 

0.171 

0.085 

482.0 

0.118 

0.059 

a" 

2 

657.6 

0.138 

0.069 

674.2 

0.093 

0.047 

FSCCSD+T*(3) 

a' 

3 

538.0 

1.391 

0.547 

550.5 

1.076 

0.431 

a' 

3 

665.9 

2.465 

0.853 

678.1 

2.872 

1.033 

a' 

2 

749.8 

2.757 

1.379 

758.4 

3.080 

1.540 

a" 

2 

459.7 

0.073 

0.037 

474.3 

0.033 

0.017 

a" 

2 

646.9 

0.044 

0.022 

662.5 

0.012 

0.006 

Given  the  data  above,  a  reasonable  estimate  of  the  error  bars  due  to  changes  in  the 
active  space  is  ±2  mH  for  both  methods.  (It  must  be  pointed  out  that  there  are  certainly 
to  be  systems  with  larger  variations,  but  more  experience  will  be  required  to  understand 
the  origin  of  the  large  variations  and  therefore  be  able  to  identify  the  worst  offenders.) 
These  correspond  to  about  ±0.05  eV  or  ±1.3  kcal  mol-1.  Although  this  may  seem  like 
a  fairly  large  uncertainty,  particularly  for  comparing  energies  at  different  points  on  a 
potential  energy  surface,  there  are  several  important  factors  worth  noting.  First,  this 
is  a  conservative  estimate  for  entirely  unknown  molecules.  Many  of  the  molecules  in 
Table  5.1  are  more  than  an  order  of  magnitude  better.    Also,  the  table  presents  only 
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the  largest  variation  —  many  others  were  significantly  smaller.  Secondly,  it  must  be 
remembered  that  for  molecules  with  symmetry,  many  of  the  IPs  exhibit  no  variation 
because  there  is  only  one  possible  choice  of  active  space  in  that  symmetry.  Finally, 
it  is  a  fairly  simple  matter,  and  not  computationally  expensive,  to  obtain  a  feel  for  the 
variation  of  any  molecule  of  interest.  Only  one  (0,0)  sector  calculation  is  required  to 
try  out  several  FSCC  active  spaces,  and  the  FSCC  calculations  themselves  are  very 
inexpensive  compared  to  the  (0,0)  reference. 

As  a  result,  while  the  uncertainties  associated  changing  the  active  space  in  T(3)  and 
T*(3)  calculations  for  the  general  case  might  be  considered  large  for  some  applications 
much  tighter  limits  may  often  be  obtained  with  minimal  effort.  For  the  remainder  of  this 
chapter,  the  active  spaces  used  will  correspond  directly  to  those  results  which  are  reported. 
Uncertainties  due  to  the  active  space  will  be  introduced  into  discussions  as  appropriate. 

5.3  Ionization  Potentials 

Although  it  can  be  used  in  many  other  ways,  the  "natural"  result  of  a  (0,1)  sector 
FSCC  calculation  is  the  energy  of  the  vertical  ionization  process  —  extracting  an  electron 
from  a  molecule  without  allowing  the  geometry  to  relax  to  that  of  the  resulting  ion.  Cal- 
culation of  ionization  potentials  is,  thus  far,  the  predominant  application  of  the  FSCCSD 
method  and  such  calculations  will  provide  the  bulk  of  the  data  for  characterization  of  the 
approximate  triples  methods.  These  results  will  be  compared  with  SRCC  calculations, 
experimental  results,  and  a  benchmark  full  CI  calculation  for  the  methylene  molecule  by 
Bauschlicher  and  Taylor  [62]. 
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Table  5.3.  Vertical  ionization  potentials  of  'A!  and  3B!  CH2  (in  eV)  calculated  with  a 
DZP  basis. 


Method 

12A, 

1A1- 

22A! 

2B2 

3B, 

2A, 

2B, 

FSCCSD 

10.21 

22.54 

14.92 

10.14 

11.30 

FSCCSD+T(3) 

10.32 

22.74 

15.00 

10.26 

11.43 

FSCCSD+T*(3) 

10.28 

22.68 

14.96 

10.23 

11.39 

Full  CI 

10.26 

22.14 

14.85 

10.16 

11.31 

QRHF-CCSD 

10.20 

22.25 

14.84 

10.14 

11.30 

QRHF-CCSD[T] 

10.25 

22.28 

14.85 

10.16 

11.31 

QRHF-CCSDT-1 

10.25 

22.28 

14.85 

QRHF-CCSDT 

10.25 

22.15 

14.85 

Note:  Geometry,  basis  set,  and  full  CI  results  taken  from  Bauschlicher  and  Taylor 
[62].  QRHF-CC  results  for  'Aj  ionizations  from  Watts,  Gauss,  and  Bartlett  [58].  3Bi 
calculations  used  a  UHF-CCSD  reference  function. 


Table  5.3  presents  FSCC  calculations  on  ionization  potentials  of  methylene  in  com- 
parison with  both  the  full  CI  data  from  Bauschlicher  and  Taylor,  and  several  QRHF-CC 
calculations.  At  the  FSCCSD  level,  four  of  the  five  calculated  IPs  are  within  0.1  eV 
of  the  full  CI  number,  with  the  'Ai-^A,  result  0.4  eV  off.  The  effect  of  the  T(3) 
correction  to  FSCCSD  is  to  move  all  of  the  results  away  from  full  CI,  but  four  of  the 
roots  are  still  within  0.1  eV,  while  the  'Aj— >22Aj  number  has  changed  to  0.6  eV  away 
from  full  CI.  Introducing  some  higher  order  effects  via  the  T*(3)  model  shifts  all  of  the 
results  back  toward  full  CI,  but  not  as  close  as  the  FSCCSD  result.  These  differences 
from  full  CI  cannot  result  from  uncertainties  due  to  the  choice  of  active  space,  which 
from  Table  5.1  are  about  0.001  eV  (0.04  mH)  for  this  particular  molecule.  By  contrast, 
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the  single  reference  QRHF-CC  calculations  start  out  with  all  five  ionizations  within  0.11 
eV  of  full  CI,  and  the  inclusion  of  triple  excitations  at  various  levels  improves  the  results 
to  match  full  CI  within  0.01  eV.  This  behavior  suggests  that  the  single  reference  CCSD 
provides  a  good  approximation  to  the  various  states  in  question,  and  that  triple  excitations 
are  relatively  unimportant,  which  would  seem  to  conflict  with  the  FSCC  results. 

It  is  important  to  recall,  however,  that  there  is  a  fundamental  difference  between 
FSCC  and  SRCC  with  respect  to  the  nature  of  their  excitation  operators.  As  discussed  in 
the  previous  chapter,  spectator  triple  excitations  must  be  added  to  the  FSCCSD  method 
to  achieve  an  expansion  space  that  is  in  some  sense  equivalent  to  that  of  the  SRCC.  But 
even  then,  the  FSCC  amplitudes  are  restricted  by  the  necessity  to  represent  several  ion- 
ized states  simultaneously.  Thus,  we  should  not  expect  the  FSCCSD  to  produce  results 
identical  to  a  single  reference  CCSD,  nor  should  we  expect  triple  excitation  corrections 
to  behave  in  the  same  way.  The  differences  between  FSCC  and  SRCC  shown  in  Table 
5.3  and  others,  below,  should  be  interpreted  as  a  manifestation  of  these  differences. 

These  differences  between  single  reference  and  Fock-space  methods  can  also  be  seen 
in  the  other  calculations  presented  here  (Tables  5.4-5.9).  The  sign  of  the  FSCC  triples 
correction  is  frequently  opposite  that  of  the  SRCC  triples  corrections,  and  their  magni- 
tudes do  not  appear  to  be  related  —  a  large  SRCC  triples  contribution  does  not  necessarily 
imply  a  large  FSCC  triples  correction,  and  vice  versa.  The  size  of  the  SRCC  and  FSCC 
triples  contributions  are  rather  different  as  well.  The  largest  SRCC  effect  is  about  0.2 
eV,  while  the  FSCCSD+T(3)  correction  is  as  large  as  0.8  eV,  and  0.5  eV  for  T*(3).  All 
this  is  further  support  for  the  idea  that  the  nature  of  the  excitation  operator  hierarchy  is 
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Table  5.4.  Vertical  ionization  potentials  of  CH2NH  (in  eV)  calculated  at  the  experimental 
geometry  with  a  DZP  basis. 


Basis/Method 

5a' 

la" 

Orbital 

4a' 

3a' 

FSCCSD 

10.34 

12.29 

15.08 

17.19 

FSCCSD+T(3) 

10.71 

12.35 

15.30 

17.46 

FSCCSD+T*(3) 

10.55 

12.29 

15.20 

17.34 

QRHF-CCSD 

10.40 

12.14 

15.17 

17.32 

QRHF-CCSDT-1 

10.43 

12.28 

15.11 

17.14 

Experiment  [63] 

10.52 

12.43 

15.13 

17.04 

Experiment  [64] 

10.70 

12.48 

15.11 

17.07 

Experiment  [65] 

10.56 

12.44 

15.00 

17.00 

Notes:  Geometries,  basis  sets,  and  QRHF-CC  results  taken  from  [66]. 

Table  5.5.  Vertical  ionization  potentials  of  CH2PH  (in  eV)  calculated  at  the  experimental 
geometry  with  a  DZP  basis. 


Basis/Method 

la" 

5a' 

Orbital 

4a' 

3a' 

FSCCSD 

10.08 

10.28 

13.17 

15.14 

FSCCSD+T(3) 

10.01 

10.33 

13.25 

15.41 

FSCCSD+T*(3) 

9.98 

10.28 

13.25 

15.27 

QRHF-CCSD 

9.90 

10.22 

13.15 

15.15 

QRHF-CCSDT-1 

10.05 

10.24 

13.15 

Experiment  [67] 

10.30 

10.70 

13.20 

15.00 

Notes:  Geometries,  basis  sets,  and  QRHF-CC  results  taken  from  [66]. 
different  in  single  reference  and  Fock-space  coupled  cluster  methods. 


75 
Table  5.6.  Vertical  ionization  potentials  for  3Sg"  O2  (in  eV)  at  the  experimental  geometry. 


Method 

Basis 

02+  Final  State 

2ng 

4V   - 

X 

4r  - 

FSCCSD 

DZP 

11.76 

17.75 

16.40 

24.16 

FSCCSD+T(3) 

DZP 

12.13 

17.99 

16.30 

24.98 

FSCCSD+T*(3) 

DZP 

12.08 

17.87 

16.29 

24.66 

UHF-CCSD 

DZP 

12.13 

17.80 

16.15 

24.64 

UHF-CCSD(T) 

DZP 

11.95 

17.77 

16.26 

24.40 

FSCCSD 

PVQZ++ 

12.37 

18.34 

16.94 

24.74 

FSCCSD+T(3) 

PVQZ++ 

12.82 

18.48 

16.74 

25.43 

FSCCSD+T*(3) 

PVQZ++ 

12.57 

18.34 

16.71 

25.11 

UHF-CCSD(T) 

PVQZ++ 

12.39 

18.26 

16.76 

24.82 

Expt.  (Adiabatic) 

12.07 

18.17 

16.11 

24.58 

Expt.  (Vertical) 

12.35 

18.33 

16.85 

24.66 

Notes:  Geometry  taken  from  Herzberg  [68].  Basis  sets,  SRCC,  and  experimental  result 
taken  from  [69]. 


Examining  the  various  tables,  we  see  that  the  SRCC  triples  correction  usually 
improves  agreement  with  experiment,  while  the  FSCC  triples  corrections  appear  to  shift 
the  calculated  result  away  from  the  experimental  one  as  often  as  towards  it.  It  is  also 
worth  noting  that  the  T*(3)  correction  is  generally  smaller  than  T(3)  and  gives  better 
agreement  with  experimental  results.  Overall,  FSCCSD  and  FSCCSD+T*(3)  seem  to 
do  about  the  same  in  comparison  with  experiment,  with  an  average  absolute  difference 
from  the  experimental  results  of  roughly  0.2  eV.  The  single  reference  CCSD  results 
have  an  average  absolute  difference  of  closer  to  0.3  eV,  which  is  improved  to  0.2  eV 
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Table  5.7.  Vertical  ionization  potentials  of  H2CO  (in  eV)  evaluated  at  the  experimental 
geometry  with  a  [5s3pldJ2slp]  basis. 


Method 


Orbital  of  ionized  electron 
2b2  lbj  5ai 


FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 

QRHF-CCSD 
QRHF-CCSD[T] 
QRHF-CCSDT-lb 
QRHF-CCSDT-3 

Experiment  [70] 


10.51 

14.36 

15.75 

10.96 

14.49 

16.30 

10.77 

14.40 

16.06 

10.59 

14.19 

15.77 

10.63 

14.37 

15.83 

10.63 

14.39 

15.83 

10.58 

14.31 

15.76 

10.88 

14.38 

16.00 

Notes:  Geometry  taken  from  Herzberg  [68].  The  basis  set  consisted  of  Dunning's 
(9s5p)— >[5.s3p]  set  on  C  and  O,  and  his  (4s)— >[2.s]  set  with  exponents  scaled  by  1.44 
on  H  [71],  with  polarization  functions  exponents  of  0.654  (C),  1.211  (O),  and  0.7  (H) 

[72]. 

by  the  addition  of  triple  excitation  effects.  Thus,  it  would  appear  that  FSCCSD  and 
FSCCSD+T*(3)  perform  about  as  well  as  single  reference  CCSD(T)  in  comparison  with 
experiment,  but  the  single  reference  calculations  (so  far)  have  an  edge  in  the  predictability. 

In  order  to  give  proper  weight  to  the  performance  of  these  methods  relative  to 
experiment,  we  must  consider  the  potential  problems  inherent  in  experimental  data,  as 
mentioned  in  the  Introduction.  Most  of  these  problems  are  displayed  by  one  or  more  of 
the  experimental  results  in  this  chapter. 

The  natural  result  of  a  FSCC  calculation  is  a  vertical  ionization  potential,  while 
experiment  may  provide  vertical,  adiabatic,  or  both.  Adiabatic  IPs  are  not  direcdy 
comparable  with  calculated  vertical  ionization  potentials,  since  they  depend  on  the  change 
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Table  5.8.  Vertical  ionization  potentials  of  N2  (in  eV)  at  /?=2.0693  au  with  a  [5s4pld] 
basis. 

Orbital  of  ionized  electron 
Method 

3<7g  1 7TU  2<7U 

FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 


15.44 

17.14 

18.64 

15.58 

16.86 

19.04 

15.44 

16.88 

18.81 

15.39 

16.69 

18.71 

15.36 

16.87 

18.57 

15.60 

16.98 

18.78 

QRHF-CCSD 

QRHF-CCSDT-1 

Experiment 

Notes:  Geometry,  basis,  and  QRHF-CC  and  experimental  results  taken  from  Rittby  and 
Bartlett  [73]. 

Table  5.9.  Vertical  ionization  potentials  of  ketene  and  diazomethane  (in  eV)  calculated 
at  experimental  geometries  calculated  with  a  DZP  basis. 


Method 


Ketene  Diazomethane 

2bi  2b2  2bi  2b2 


9.40 

14.25 

8.60 

13.79 

9.39 

14.49 

8.52 

13.85 

9.35 

14.38 

8.47 

13.77 

9.29 

14.31 

8.43 

13.78 

9.38 

14.24 

8.54 

13.73 

FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 

QRHF-CCSD 
QRHF-CCSD(T) 


Experiment 9.64  (a)        13.84  (a)         9.00  (a)        14.13  (v) 

Notes:    Geometries,  basis  sets,  and  experimental  results  taken  from  Rittby,  Pal,  and 
Bartlett  [55].  Experimental  IPs  are  adiabatic  or  vertical,  as  indicated  in  the  table. 

in  geometry  from  the  neutral  to  the  ion.  If  the  ion  and  neutral  geometries  are  very  similar, 
the  vertical  and  adiabatic  IPs  will  also  be  nearly  the  same,  but  in  general,  the  adiabatic 
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IP  is  only  a  lower  bound  for  the  vertical  one.   In  cases  where  the  experimental  value 
is  adiabatic  or  unspecified,  a  larger  difference  from  the  computed  vertical  IP  must  be 
expected. 

A  second  problem  with  experimental  results  is  that,  as  with  theoretical  methods, 
different  experimental  techniques  can  give  somewhat  different  results.  This  is  illustrated 
in  Table  5.4,  which  lists  three  sets  of  experimental  data  for  methyleneamine,  differing 
among  themselves  by  as  much  as  0.2  eV. 

A  third  complication  of  using  experimental  data  is  that  details  of  the  calculation 
other  than  the  correlation  method  become  important,  particularly  the  basis  set.  Most 
of  the  calculations  presented  here  use  a  DZP  or  slightly  better  quality  basis  set.  While 
very  useful  for  comparing  one  theoretical  method  to  another,  these  basis  sets  may  not  be 
adequate,  however,  to  obtain  good  results  compared  to  experiment.  Table  5.6  shows  the 
ionization  potentials  of  O2  as  calculated  using  both  DZP  (32  functions)  and  an  augmented 
PVQZ  basis  (126  functions)  in  order  to  demonstrate  the  effect  of  basis  set  on  both 
the  FSCC  and  SRCC  IPs.  It  is  clear  that  the  PVQZ++  basis  provides  a  much  better 
comparison  with  the  experimental  results  than  the  DZP  basis,  for  both  FSCC  and  SRCC 
calculations.  The  change  in  FSCC  results  due  to  increasing  the  basis  set  is  about  0.5-0.6 
eV,  and  for  the  SRCC  0.4-0.5  eV.  Although  there  is  no  experimental  comparison,  a 
similar  basis  set  effect  can  also  be  seen  in  the  HOF  results  in  Table  5.2,  where  the 
change  from  a  DZP  to  a  PVQZ  basis  changes  the  IPs  by  0.2-0.4  eV.  It  would  seem, 
then,  that  correlated  calculations  with  DZP-quality  basis  sets  should  not  be  expected  to 
reproduce  experimental  results  too  well. 
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Given  that  the  basis  sets  used  here  are  generally  too  small,  coupled  with  the  other 
problems  of  the  experimental  data  themselves,  it  would  seem  that  an  error  of  0.2  eV  for 
these  calculations  is  not  unreasonable,  and  is  not  necessarily  an  indication  of  intrinsic 
problems  with  the  method.  Indeed,  the  fact  that  the  SRCC  results  show  a  similar  margin 
of  error  with  respect  to  experiment  is  a  good  sign  for  FSCC.  It  would  seem  that  the 
biggest  difference  in  this  regard  between  the  single  reference  and  Fock-space  results  is 
the  relative  lack  of  consistency  of  the  FSCC  triples  corrections,  which  can  probably  be 
better  understood  with  more  calculations. 

5.4  Potential  Energy  Surfaces  and  Related  Properties 

Another  use  of  FSCC  methods  is  for  the  study  of  potential  energy  surfaces,  especially 
with  the  ability  to  look  at  several  surfaces  in  the  same  calculation.  The  basic  result  of 
the  FSCC  calculation  is  still  an  energy  difference,  but  when  added  to  the  absolute  energy 
of  the  (0,0)  sector  reference,  it  provides  an  absolute  energy  for  the  state  of  interest. 
In  the  (0,1)  sector  case,  the  state  of  interest  has  one  electron  removed  relative  to  the 
reference  (0,0)  state.  While  this  may  seem  a  round  about  way  of  calculating  the  PES 
of  a  molecule,  it  has  practical  advantages  in  many  cases.  For  calculations  on  molecules 
with  low-lying  electronic  states  or  symmetry  breaking  phenomena,  direct  calculation  of 
the  PES  by  traditional  single  reference  methods  frequently  suffers  from  convergence 
and  other  problems  due  to  the  proximity  of  other  states  to  the  one  of  interest.  On  the 
other  hand,  the  same  molecule  with  another  electron  is  often  well-behaved  at  the  same 
geometries.    This  the  same  as  the  ionization  process,  but  from  the  opposite  point  of 
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view  —  in  this  case,  we  want  to  talk  about  the  final  state,  while  in  IPs,  it  is  traditional 
to  look  at  it  from  the  point  of  view  of  the  initial  state. 

This  application  of  FSCC  methods  is  a  very  recent  development,  and  there  are  only 
a  few  such  applications  in  the  literature.  Moreover,  because  they  deal  with  symmetry 
breaking  problems,  there  is  little  conclusive  experimental  evidence  to  allow  concrete 
conclusions  to  be  drawn  about  the  performance  theoretical  methods.  Comparisons  are 
possible  with  single  reference  calculations,  though  as  we  shall  see,  the  comparison 
of  high-level  SRCC  calculations  with  Fock-space  results  does  not  always  lead  to  an 
unambiguous  determination  of  which  approach  is  more  reliable.  Nevertheless,  FSCC  is 
very  well  suited  to  such  applications  and  will  undoubtedly  see  increasing  use  in  the  study 
of  potential  energy  surfaces  (beginning  with  the  next  chapter  of  this  work),  so  it  is  worth 
trying  to  get  a  feel  for  its  performance. 

Engelbrecht  and  Liu,  in  a  1983  paper  [74],  used  multiconfiguration  self-consistent 
field  (MCSCF)  plus  CI  calculations  to  characterize  the  3A2  state  of  CO2.  This  state, 
which  is  bent  (C2V  symmetry),  has  the  SCF  occupation 


">rvl'>  ,       ^rul  r       9.,2, 


(core)3a J 2^24^1 369 5a jl6jl 024^2601 
and  interacts  strongly  with  the  configuration 

because  both  the  a,  and  bj  orbitals  transform  as  a"  in  C,  symmetry.   As  a  result,  SCF 
calculations  can  obtain  a  lower  energy  by  going  to  a  Cs  structure. 
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Asymmetric  Stretch  Potential  of  A2  CO- 


-0.10 


-0.05 


0.00 

AR(A) 


0.05 


FSCCSD 


FSCCSD+T(3) 
FSCCSD+T*(3) 
QRHF-CCSD 
QRHF-CCSD(T) 


0.10 


Figure  5.1.  Potential  curves  for  the  asymmetric  stretch  mode  (AR  =  ^(Ri  —  Ro))  of 
3A2  CO2  at  various  levels  of  theory.  The  C2V  structure  is  the  DZP/FOCIB  geometry  of 
Engelbrecht  and  Liu  [74]. 

With  suitable  methods,  such  interactions  can  be  treated  properly,  and  reason- 
able potential  energy  surfaces  obtained,  but  this  requires  the  ability  to  describe 
inherently  multireference  states.  Most  such  methods  are  not  many-body  in  na- 
ture, and  of  course  we  prefer  a  size  extensive  approach  if  possible  (Chapter 
3). 

It  has  been  observed  that  QRHF-based  coupled  cluster  methods  are  capable  of  han- 
dling many  multireference  problems  when  sufficient  triple  excitation  effects  are  incor- 
porated [34].  In  Figure  5.1,  we  see  two  QRHF-CC  potential  curves  for  the  asymmet- 
ric stretch  of  CO2,  which  were  calculated  from  closed-shell  C022~  SCF  orbitals.  The 
QRHF-CCSD  curve  has  a  double  minimum  shape,  while  the  QRHF-CCSD(T)  method 
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gives  a  stable  C2V  structure.    This  is  consistent  with  the  MCSCF+CI  calculations  of 
Engelbrecht  and  Liu  and  their  conclusion  that  the  equilibrium  structure  of  this  state  is 

C2v 

The  results  of  three  sets  of  Fock-space  coupled  cluster  calculations  using  2A\  CC>2~ 
as  the  (0,0)  reference  are  also  displayed  in  Figure  5.1.  FSCCSD  gets  the  proper  shape 
of  the  PES  even  without  triple  excitation  effects,  and  the  inclusion  of  triples  makes 
small  adjustments  to  the  shape  of  the  surface.  (Note  that,  as  in  the  IP  calculations,  the 
FSCCSD+T*(3)  results  lie  between  FSCCSD  and  FSCCSD+T(3).)  This  behavior  is  what 
we  might  expect  from  a  method  that  is  designed  in  a  multireference  framework. 

Another  symmetry-breaking  problem  to  which  a  number  of  single  and  multiref- 
erence approaches  have  been  applied  is  the  nitrate  radical,  NO3.  Different  exper- 
iments have  been  interpreted  in  terms  of  either  a  totally  symmetry  D3h  equilib- 
rium geometry,  or  a  lower  symmetry  C2V  structure  with  one  N-0  bond  length  dif- 
ferent from  the  other  two.  However,  as  Kawaguchi  and  coworkers  [75]  comment, 
no  one  has  advanced  a  model  which  explains  all  of  the  available  spectroscopic 
data. 

In  the  most  recent  theoretical  examination  of  the  problem,  Stanton,  Gauss,  and  Bartlett 
[76]  located  three  stationary  points  on  the  NO3  surface  using  the  QRHF-CCSD  method. 
The  D3h  structure  was  determined  to  be  a  local  minimum,  and  two  C2v  structures  were 
found.  They  are  characterized  by  the  relationship  between  the  unique  N-0  bond  distance 
and  the  other  two:  one  long/two  short  (1L2S)  or  one  two  long/one  short  (2L1S).  The  1L2S 
structure  was  also  found  to  be  a  local  minimum,  and  the  2L1S  a  transition  state  192  cm"1 
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Table  5.10.  Calculated  relative  energies  (in  kcal  mol"1)  of  Civ  and  D3h  NO3  structures. 


Method 

DZP  basis 
1L2S-D3h            2L1S-D3h 

PVTZ  basis 
1L2S-D3h            2L1S-D3h 

FSCCSD 

2.50 

1.05 

2.94 

1.05 

FSCCSD+T(3) 

3.03 

1.64 

3.72 

1.73 

FSCCSD+T*(3) 

2.74 

1.35 

3.36 

1.41 

QRHF-CCSD 

-2.59 

-2.04 

-2.96 

QRHF-CCSD(T) 

0.34 

0.57 

QRHF-CCSDT-3 

0.21 

Notes:  Geometries,  optimized  at  the  DZP/QRHF-CCSD  level,  taken  from  [76].  The  C2v 
geometries  are  characterized  by  one  of  the  NO  bonds  being  longer  (1L2S)  or  shorter 
(2L1S)  than  the  other  two.  QRHF-CC  relative  energies  taken  from  [77]. 

higher  in  energy  than  the  1L2S  state;  both  C2V  structures  being  lower  in  energy  than  the 
totally  symmetry  structure.  Their  conclusion  is  that  the  molecule  has  a  C2V  equilibrium 
geometry,  but  probably  has  an  extremely  low  barrier  to  pseudorotation,  which  would  give 
spectra  with  characteristic  of  a  D3h  system. 

In  a  later  paper  focusing  on  different  orbital  choices  for  single  reference  treatments  of 
symmetry  breaking  problems,  Stanton  e t  al.  [77]  calculated  relative  energies  of  the  D3h 
and  C2v  minima  using  more  sophisticated  methods  than  in  their  first  paper.  Their  QRHF- 
CC  results  are  shown  in  Table  5.10  along  with  FSCC  results  using  the  same  basis  sets 
and  geometries.  Using  both  DZP  and  PVTZ  basis  sets,  QRHF-CCSD  calculations  favor 
the  C2v  minimum  by  about  3  kcal  mol"1.  As  soon  as  triples  corrections  are  introduced, 
either  using  the  noniterative  CCSD(T)  approximation  or  the  iterative  CCSDT-3  model, 
the  energies  shift  in  favor  of  the  totally  symmetric  structure,  but  by  only  0.5  kcal  mol"1. 
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FSCC,  on  the  other  hand,  predicts  the  D3h  structure  to  be  lower  for  all  three  methods 
shown,  but  by  roughly  3  kcal  mol"1. 

Given  previous  experience  with  QRHF-CC  methods  applied  to  problems  of  this  sort, 
we  should  have  greater  confidence  in  the  QRHF-CCSD(T)  and  T-3  results  than  in  results 
without  triple  excitations,  but  with  the  available  data,  it  is  impossible  to  tell  whether 
the  single  reference  numbers  including  triple  excitations  or  the  FSCC  results  are  more 
reliable.  As  Stanton  et  al.  [77]  observe,  "An  extremely  sophisticated  (and  expensive) 
calculation  would  be  required  to  firmly  establish  the  nature  of  the  D3h  stationary  point." 

C3+  is  another  molecule  for  which  both  experimental  and  theoretical  studies  have 
produced  differing  opinions  on  the  nature  of  the  ground  state.  Coulomb  explosion  data  can 
be  interpreted  in  terms  of  either  a  floppy  linear  molecule  or  a  bent  structure,  depending  on 
the  vibrational  temperature  of  the  sample  [78,79].  On  the  theoretical  side,  the  consensus 
is  that  the  ground  state  is  a  bent  2B2  with  the  linear  2EU+  higher  in  energy  [80].  The 
energy  difference  between  these  two  states  is  still  not  a  settled  issue.  This  is  important 
because  if  the  barrier  is  small  enough,  the  molecule  may  be  quasilinear  —  although  the 
true  ground  state  is  bent,  a  large  amplitude  bending  motion  could  bring  the  molecule 
through  the  linear  structure. 

Tables  5.1 1  and  5.12  present  the  results  of  FSCC  optimizations,  vibrational  frequency 
calculations,  and  relative  energies  in  comparison  with  single  reference  CC  results  from 
Watts,  Stanton,  Gauss,  and  Bartlett  [80]  which  include  full  CCSDT  calculations.  All  three 
Fock-space  methods  give  essentially  the  same  geometry  for  both  structures,  comparing 
most  closely  with  the  single  reference  CCSD  results.  For  the  vibrational  frequencies  of  the 
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Table  5.11.  Optimized  geometries  and  vibrational  frequencies  of  2B2  C3+,  using  a  DZP 
basis. 


Method 

R(A) 

on 

u>i  (ai) 
(cm"1) 

lj2  (aj) 
(cm"') 

u>3  (b2) 
(cm"') 

FSCCSD 

1.328 

71.6 

1685 

583 

1269 

FSCCSD+T(3) 

1.336 

70.2 

1643 

614 

1221 

FSCCSD+T*(3) 

1.335 

70.5 

1651 

603 

1225 

CCSD 

1.337 

68.3 

1668 

687 

1215 

CCSD(T) 

1.350 

68.3 

1601 

515 

1074 

CCSDT-1 

1.350 

69.0 

1558 

569 

1078 

CCSDT-2 

1.347 

69.0 

1612 

595 

1213 

CCSDT 

1.350 

68.0 

1601 

638 

1021 

Notes:  Basis  set  and  SRCC  results  taken  from  Watts,  et  al.   [80]. 

Table  5.12.  Optimized  geometries,  and  vibrational  frequencies  of  2EU+  C3+,  and  energy 
relative  to  2B2  C3+,  using  a  DZP  basis. 


Method 

R(A) 

^(7TU) 

(cm1) 

^'(<7g) 

(cm"') 

(cm"') 

E(2SU+  - 

2B2) 

(kcal  mol') 

FSCCSD 

1.312 

88 

1227 

492/ 

7.1 

FSCCSD+T(3) 

1.315 

71/ 

1203 

1911 

10.9 

FSCCSD+T*(3) 

1.316 

532 

1205 

1462 

8.6 

CCSD 

1.314 

2500/ 

13.2 

CCSD(T) 

1.327 

1431/ 

2.9 

CCSDT-1 

1.339 

1135 

-1.7 

CCSDT-2 

1.329 

1595 

2.2 

CCSDT 

1.327 

451/ 

4.9 

Notes:  Basis  set  and  SRCC  results  taken  from  Watts,  et  al.   [80]. 

2B2  state,  the  FSCC  results  are  also  basically  the  same,  and  correspond  most  closely  with 
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the  CCSDT-2  approximate  iterative  triples  model  among  the  single  reference  calculations. 
On  the  basis  of  several  calculations  with  a  larger  basis  set,  Watts  and  coworkers  suggest 
that  the  true  harmonic  frequencies  should  lie  in  the  ranges  1590-1629,  592-638,  and 
1094-1173  cm"1,  and  compared  to  these,  the  FSCC  u\  and  u>3  are  slightly  high. 

The  asymmetric  stretching  frequency  (<ru)  of  the  2EU+  state  (Table  5.12)  is  a  more 
complex  problem.  An  imaginary  frequency  in  this  mode  indicates  that  the  D^  structure 
is  a  transition  state,  and  a  lower  energy  can  be  obtained  by  changing  the  bond  lengths 
asymmetrically.  For  both  single  reference  and  Fock-space  methods,  the  calculated 
frequency  varies  widely.  Also,  the  bending  mode  (iru)  changes  dramatically  with  the 
choice  of  FSCC  method.  Such  large  changes  due  to  different  correlation  methods  are 
usually  an  indication  that  important  aspects  of  the  problem  are  not  being  taken  into 
account,  and  a  more  detailed  examination,  using  both  single  reference  and  Fock-space 
methods  would  be  in  order. 

Table  5.12  also  presents  the  energy  differences  between  the  linear  and  bent  structures. 
Note  that  the  single  reference  CCSDT-1  method  places  the  linear  structure  lower  in 
energy  than  the  bent  one.  Watts  et  al.  attributed  this  to  an  apparent  tendency  for  the 
iterative  nature  of  the  CCSDT-1  method  to  overemphasize  problems  in  the  closely-related 
CCSD(T)  model  in  difficult  cases,  thus  producing  poorer  results  than  the  perturbative 
CCSD(T)  method. 

Their  CCSDT  AE  of  4.9  kcal  mol"1  is  in  reasonable  accord  with  CISDTQ  calculations 
of  Grev,  Alberts,  and  Schaefer  [81,82]  indicating  a  4.1  kcal  mol"1  difference.  Taylor  et 
al.   [83]  used  multireference  CI  methods  to  obtain  a  separation  of  1.68  kcal  mol"1  with 
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a  DZP  basis  set,  and  claim  that  this  should  be  very  close  to  the  full  CI  limit.  These 
authors  also  show  that  multireference  effects  are  important  for  both  the  bent  and  linear 
structures.  Despite  the  multireference  character  of  the  FSCC  approach,  all  three  methods 
presented  here  give  a  significantly  larger  energy  gap  than  Taylor's  calculations.  The 
problem  with  the  relative  energies  of  the  two  states  may  be  due,  at  least  in  part,  to  the 
apparent  problems  in  properly  describing  the  2SU+  state.  It  would  seem,  on  this  minimal 
evidence,  that  where  single  reference  methods  including  triple  excitations  have  problems, 
Fock-space  methods  do  as  well. 

5.5  Conclusions 

A  variety  of  results  has  been  presented  using  the  Fock-space  coupled  cluster  methods 
FSCCSD,  FSCCSD+T(3),  and  FSCCSD+T*(3),  for  both  ionization  potential  problems  at 
fixed  geometries,  and  in  cases  where  the  potential  surfaces  must  be  explored. 

Overall,  the  performance  of  the  FSCC  methods  is  good.  FSCCSD  results  usually 
compare  well  with  single  reference  CCSD(T)  calculations.  The  FSCC  triples  corrections 
seem  initially  to  be  less  predictable,  than  the  triples  correction  in  single  reference  methods. 
The  Fock-space  triple  excitations  have  a  different  nature  than  their  SRCC  counterparts, 
and  it  will  require  more  experience  with  these  methods  to  understand  them  and  discern 
the  system  they  follow. 

Including  some  higher-order  contributions  in  the  T*(3)  method  seems  to  have  a 
"dampening"  effect  compared  to  the  basic  T(3)  model,  which  causes  the  T*(3)  results  to 
compare  better  with  experimental  results  or  other  calculations.  This  is  probably  similar 
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to  the  single  reference  coupled  cluster  case,  where  the  simplest  perturbative  inclusion 
of  triples,  CCSD+T(CCSD)  was  found  to  be  significantly  improved  by  the  inclusion 
of  certain  higher  order  terms,  resulting  in  the  CCSD(T)  model.  For  the  present,  the 
FSCCSD+T*(3)  method  seems  to  be  preferable  to  the  purely  third-order  FSCCSD+T(3). 
In  the  long  run,  it  is  probably  worthwhile  to  investigate  more  systematic  ways  of 
incorporating  higher-order  terms  in  the  FSCC  triples  contribution,  as  discussed  in  the 
previous  chapter. 

Numerically,  effects  of  changing  the  active  space  in  the  FSCCSD+triples  methods 
have  been  demonstrated  to  introduce  uncertainties  as  large  as  ±2  mH  in  calculated 
energies  (IPs  or  state  energies),  however  most  of  the  molecules  examined  have  much 
lower  uncertainties.  On  the  whole,  both  FSCCSD  and  FSCCSD+T*(3)  seem  to  do  about 
as  well  compared  to  experiment  as  SRCC  methods  including  perturbative  triples,  and 
better  than  CCSD.  For  very  accurate  comparison  with  experiment  for  both  FSCC  and 
single  reference  methods,  however,  a  large  basis  is  required. 

For  energy  differences  on  potential  energy  surfaces,  where  the  requirements  are  tighter 
than  on  IPs,  Fock-space  results  do  not  seem  to  compare  as  well  with  single  reference 
methods  as  for  IPs.  Although  the  signs  of  relative  energies  evaluated  using  FSCC  seem 
to  match  higher-quality  single  reference  calculations,  the  Fock-space  energy  gaps  are 
noticeably  larger  than  in  the  SRCC  case.  The  cases  examined  are  hard  problems,  both 
experimentally  and  theoretically,  and  it  is  not  clear  that  either  approach  can  be  claimed 
to  be  better.  Consequently,  great  care  must  be  taken  with  such  results,  and  more  study 
would  be  most  useful. 
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The  FSCC  methods  all  seem  to  do  quite  well  in  determining  local  properties  of 
potential  surfaces  —  geometries,  vibrational  frequencies,  etc.  In  most  of  these  cases, 
there  does  not  appear  to  be  a  substantial  difference  among  the  three  FSCC  methods 
examined,  except  for  linear  C3+,  for  which  both  single  reference  methods  including  triple 
excitations  and  Fock-space  methods  seem  to  have  problems. 

In  addition  to  the  quality  of  results  of  a  new  method,  one  should  also  consider  its  ease 
of  use  and  cost  in  comparison  with  competitive  methods.  By  this  measure,  FSCC  methods 
have  a  decided  advantage  over  the  single  reference  approaches.  An  FSCC  calculation, 
even  including  triple  excitation  effects,  amounts  to  a  small  fraction  of  the  cost  of  the 
(0,0)  sector  CCSD  reference  calculation.  In  contrast,  each  of  the  IP  calculations  required 
separate  SRCC  calculations  for  the  neutral,  and  for  each  ion  state,  and  at  least  one  of 
those  on  an  open-shell  molecule,  for  which  the  SRCC  is  about  four  times  more  expensive 
than  a  closed  shell  calculation.  In  most  cases,  the  reference  CCSD  calculation  used  in 
the  Fock-space  approach  is  a  closed  shell,  and  all  of  the  IPs  can  be  computed  in  a  single 
calculation.  For  potential  energy  surfaces,  the  ability  to  explore  several  surfaces  at  once 
is  also  a  distinct  advantage  over  single  reference  methods. 


CHAPTER  6 
THE  STRUCTURE  OF  THE  C5+  MOLECULE 


6.1   Background 

Although  there  is  relatively  little  experimental  data  on  C3+,  discussed  in  the 
preceding  chapter,  there  is  even  less  known  about  Cs+.  The  experiments  which 
have  observed  Cs+  have  done  so  with  a  mass  spectrometer  in  reactivity  [84],  frag- 
mentation [85-87],  or  mobility  [88]  studies.  The  essential  conclusion  of  these  pa- 
pers (with  respect  to  Cs+)  is  that  there  is  a  single  isomer  (multiple  isomers  show 
up  starting  with  C7+)  and  that  its  reactivity  and  mobility  are  similar  to  those  of 
its  neighbors,  C4+  and  C6+,  and  the  whole  set  of  clusters  from  C2+  to  C6+  are 
commonly  labeled  as  "linear"  molecules  [84,88].  (C3+,  which  theoretical  results 
agree  is  actually  bent,  fits  the  trends  for  "linear"  clusters  in  the  experimental  data 
as  well  as  Cs+  does,  suggesting  that  the  "linear"  label  cannot  be  interpreted  rigor- 
ously.) 

Recently,  matrix-isolated  infrared  spectroscopy  experiments,  using  isotopically 
mixed  graphite,  produced  an  IR  band  that  could  be  attributed  to  a  planar  pentago- 
nal (Dsh  symmetry)  Cs+  molecule  [89].  The  neutral  C5  cluster  has  been  well  es- 
tablished, both  theoretically  and  experimentally,  to  have  a  linear  ground  state  [90], 
with  cyclic  structures  much  higher  in  energy.  In  this  situation,  it  would  seem  un- 
likely that  the  cation  would  be  so  different  from  the  neutral  that  a  cyclic  cation 
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would  be  observable,  but  the  existence  of  isoenergetic  linear  and  cyclic  structures 
of  neutral  C4,  which  is  now  generally  accepted  [90,91],  was  once  also  thought  un- 
likely. 

A  series  of  SCF  and  MBPT(2)  calculations  [92]  established  that  Dsh  Cs+  could  not 
be  responsible  for  the  observed  band,  and  a  reanalysis  of  the  experimental  data  led  to 
assignment  of  the  band  to  a  C3»H20  adduct,  in  agreement  with  a  later  assignment  by 
Ortman  et  al.  [93].  The  low-level  calculations  pointed  toward  a  linear  ground  state 
for  C5+,  but  there  were  major  problems  with  spin  contamination  of  the  UHF-based 
calculations,  stability  of  the  SCF  solution,  symmetry  breaking,  and  failure  to  converge 
calculations  on  some  states.  In  this  chapter  we  attempt  to  understand  the  source  of  the 
troubles  with  lower-level  calculations  and  to  establish  the  true  ground  state  structure  of 
C5+  using  high-quality  many-body  methods. 

6.2  Computational  Methods 

The  majority  of  the  calculations  in  this  work  used  a  DZP  basis  set  comprised  of 
Dunning' s  (9s5p)/[4s2p]  double-zeta  set  [71]  and  a  polarization  function  (a=0.654)  as 
optimized  by  Redmon,  Purvis,  and  Bartlett  [72].  All  six  cartesian  d  components  were 
retained,  giving  a  total  of  80  contracted  basis  functions.  At  selected  points  on  the 
various  potential  energy  surfaces,  a  larger  PVTZ  basis  [61]  was  used  to  compute  the 
energy.  This  is  a  {\Qs5p2d\f)l[As?>p2d\f\  generally-contracted  basis,  which  was  used  with 
spherical  harmonic  polarization  components  (5 J  and  If),  for  a  total  of  150  contracted 
basis  functions. 
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Both  single  reference  QRHF-based  and  Fock-space  coupled  cluster  methods  have 
been  used  and  in  both  types  of  calculations,  triple  excitation  effects  have  been  con- 
sidered, mostly  through  perturbative  corrections:  QRHF-CCSD(T),  FSCCSD+T(3),  and 
FSCCSD+T*(3).  Both  QRHF-CC  and  FSCC  calculations  used  the  closed  shell  neutral  C5 
molecule  as  a  reference.  The  reasons  for  using  both  SRCC  and  Fock-space  approaches 
are  both  historical  and  practical.  When  work  on  this  project  was  first  begun,  no  effi- 
cient FSCC  program  was  available  and  QRHF-based  approaches  were  the  best  choice  to 
deal  with  the  problems  observed  in  the  UHF  calculations  mentioned  earlier.  When  the 
newly-completed  ACES  II  FSCCSD  code  was  applied  to  this  molecule,  the  nature  of  the 
problems  became  apparent,  and  it  was  clear  that  understanding  the  Cs+  molecule  was 
a  problem  well-suited  to  this  method  because  of  its  closely-spaced,  strongly  interacting 
electronic  states. 

In  some  calculations,  the  carbon  \s  orbitals  and  the  corresponding  virtual  orbitals  (five 
occupied  and  five  virtual  for  the  DZP  basis  with  cartesian  angular  momentum  functions, 
five  occupied  and  four  virtual  for  the  spherical  harmonic  PVTZ  basis)  were  not  correlated 
in  order  to  reduce  the  computational  cost.  The  effect  of  dropping  these  orbitals  on  the 
relative  energies  is  negligible. 

Geometry  optimizations  have  been  carried  out  using  QRHF-CCSD  and  all  three 
Fock-space  methods.  Analytic  gradient  techniques  were  used  with  the  QRHF-CCSD 
calculations,  but  gradients  have  not  been  formulated  for  FSCC  yet,  so  optimizations 
were  performed  via  numerical  differentiation  of  FSCC  energies.  Similarly,  vibrational 
frequencies  were  evaluated  from  the  numerical  second  derivatives  of  the  FSCC  energy, 
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or  by  numerical  differentiation  of  the  QRHF-CCSD  analytic  gradients.  All  optimizations 
and  frequency  calculations  used  the  DZP  basis  set. 

6.3  The  Heart  of  the  Problem 

In  neutral  C5,  the  highest  occupied  orbital  in  the  Hartree-Fock  approximation  is 
of  7Tg  symmetry,  and  the  next  highest  are  au+  and  og~.  Using  Koopmans'  theorem 
ionization  potentials,  the  two  2S  states  of  Cs+  obtained  by  removing  electrons  from 
either  of  the  a  orbitals  should  be  more  than  60  kcal  mol"1  above  the  2n.  Unfortunately, 
for  molecules  not  well-described  at  the  SCF  level,  as  the  original  calculations  suggested 
for  C5+,  Koopmans'  approximation  is  not  very  useful. 

QRHF-CCSD  optimization  of  the  2ng  state  went  well,  but  there  were  problems  with 
the  S  states:  The  2EU+  state  optimized  to  a  stationary  point,  but  only  after  significant 
computational  effort  because  large  T,  and  T2  amplitudes  made  the  calculations  very 
hard  to  converge,  and  the  "Eg"  state  could  not  be  converged  at  all.  In  addition, 
attempts  to  calculate  the  vibrational  frequencies  of  the  2EU+  state,  to  see  if  it  is  in 
fact  a  local  minimum,  produced  unphysical  results  for  the  au  modes.  Such  problems 
are  characteristic  of  a  symmetry  breaking  or  state  interaction  problem  —  the  large 
cluster  amplitudes  indicating  that  other  determinants  mix  strongly  with  the  SCF  reference 
determinant,  and  the  unphysical  vibrational  frequencies  a  reflection  of  a  discontinuity 
in  the  PES,  often  caused  by  individual  calculations  in  the  finite  difference  frequency 
evaluation  that  may  have  converged  to  a  low-lying  excited  state  relative  to  the  one  of 
interest. 
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Potential  Energy  Surfaces  of  Linear  C 
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Figure  6.1.  DZP/FSCCSD  potential  curves  for  three  states  of  Cs+,  centered  on  the 
DZP/QRHF-CCSD  2SU+  stationary  point. 

A  FSCCSD  calculation  at  the  2SU+  geometry  revealed  that  all  three  states  were  within 

about  2.2  kcal  mol"1  rather  than  the  60  kcal  mol"1  separation  predicted  by  Koopmans' 

theorem.  The  proximity  of  the  two  £  states,  which  couple  by  an  asymmetric  stretching 

(<7U)  motion  produces  a  pseudo-Jahn-Teller  effect  (PJTE),  giving  the  2SU+  PES  a  minimum 

at  a  Coov  structure,  breaking  the  perpendicular  symmetry  plane  and  the  center  of  inversion. 

This  is  shown  in  Figure  6.1,  which  displays  the  energy  of  the  three  states  for  motion  along 

the  au  mode  leading  to  the  minimum  energy  for  the  2SU+  state.  Such  a  system  is  ideally 

suited  for  the  application  of  FSCC  (0,1)  methods  since  not  only  do  the  states  require  a 

multireference  description,  but  the  neutral  C5  molecule  provides  a  well-behaved,  closed 

shell  (0,0)  sector  reference  function. 
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Table  6.1.  Optimized  geometries  of  linear  stationary  points  of  Cs+  (in  A). 


State/Method 

Ci  -  C2  -  C3  - 

C4  —  C5 

2ng  c5+ 

R12  (=R45) 

R23  (=R34) 

QRHF-CCSD 

1.351 

1.299 

FSCCSD 

1.350 

1.300 

FSCCSD+T(3) 

1.352 

1.301 

FSCCSD+T*(3) 

1.353 

1.301 

2S+  c5+ 

R12 

R23 

R34 

R45 

QRHF-CCSD 

1.245 

1.347 

1.265 

1.344 

FSCCSD 

1.237 

1.350 

1.258 

1.348 

FSCCSD+T(3) 

1.254 

1.327 

1.274 

1.330 

FSCCSD+T*(3) 

1.248 

1.339 

1.265 

1.339 

2su+  c5+ 

Ri: 

!  (=R45) 

R23  (=R34) 

QRHF-CCSD 

1.288 

1.299 

FSCCSD 

1.291 

1.299 

FSCCSD+T(3) 

1.295 

1.299 

FSCCSD+T*(3) 

1.294 

1.299 

C5a) 

R.: 

!  (=R45) 

R23  (=R34) 

105  CGTO/CEPA-1 

1.289 

1.283 

H 

—  C'i  —  C2  —  C3 

—  C4  —  C5 

HC5+b) 

Rl2 

R23 

R34 

R45 

111  CGTO/CEPA-1 

1.227 

1.335 

1.254 

1.336 

a)  From  Botschwina  and  Sebald  [94]. 
h)  From  Botschwina  [95];  Rm  =  1.076  A. 

6.4  Geometries  of  the  Linear  Stationary  Points  of  Cs+ 

Once  the  origin  of  the  convergence  and  other  problems  with  the  QRHF-CC  cal- 
culations was  recognized,  the  situation  was  greatly  improved.    The  2Sg"  state  could 
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immediately  be  eliminated  from  consideration  as  a  possible  ground  state  of  the  system, 
and  symmetry  lowering  by  the  2EU+  state  due  to  the  PJTE  interaction  also  increases 
its  separation  from  the  2ng  state,  which  should  reduce  second-order  vibronic  coupling 
between  the  two  states. 

Table  6.1  presents  the  final  optimized  geometries  of  the  two  Dod,  and  one  Coov 
linear  stationary  points  using  both  QRHF-CC  and  FSCC  methods.  For  comparison, 
the  geometries  of  neutral  C5,  as  optimized  by  Botschwina  and  Sebald  [94],  and  the 
closely  related  system,  HCs"1",  as  optimized  by  Botschwina  [95].  The  three  Fock-space 
approaches  and  QRHF-CCSD  all  agree  on  the  structures  of  these  stationary  points, 
with  differences  among  methods  similar  to  those  seen  for  C3+  in  the  previous  chap- 
ter, and  certainly  no  larger  than  would  be  expected  for  different  correlation  meth- 
ods. 

The  2£u+  structure  is  very  similar  to  the  neutral  C5  geometry,  having  nearly 
equal  bond  lengths  with  the  inner  one  slightly  longer.  The  2ng  state  has  longer 
outer  bonds,  and  a  more  pronounced  difference  between  inner  and  outer  bond  dis- 
tances. The  Cqov  2S+  structure  is  essentially  a  shift  of  atoms  2  and  4  towards 
1  and  3,  respectively,  from  the  2EU+  structure,  resulting  in  a  molecule  with  alter- 
nating bond  lengths  but  essentially  the  same  overall  length  (at  the  FSCCSD+T*(3) 
level,  it  lengthens  by  0.1  A,  but  for  the  others  the  difference  is  insignificant). 
Interestingly,  the  2S+  structure  is  also  very  similar  to  that  of  HCs+,  suggest- 
ing the  possibility  that  the  two  species  might  have  other  common  characteris- 
tics. 
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Table  6.2.  Relative  energies  (in  kcal  mol"1)  of  linear  Cs+  stationary  points. 


Method 

Structure3* 

Energy  1 

2ng 

Relative  to  2S+ 

2Y  + 

DZP/QRHF-CCSD 

QRHF-CCSD 

11.9 

29.9 

DZP/QRHF-CCSD[T] 

QRHF-CCSD 

6.1 

3.9 

DZP/FS  CCSD 

FS  CCSD 

0.6 

6.7 

DZP/FS  CCSD+T(3) 

FS  CCSD+T(3) 

-13.9 

2.0 

DZP/FS  CCSD+T*(3) 

FS  CCSD+T*(3) 

-9.1 

3.9 

PVTZ/QRHF-CCSD 

QRHF-CCSD 

13.4 

30.1 

PVTZ/QRHF-CCSD[T] 

QRHF-CCSD 

8.6 

3.8 

PVTZ/FS  CCSD 

FS  CCSD 

3.1 

6.9 

PVTZ/FS  CCSD+T(3) 

FS  CCSD+T(3) 

-11.7 

5.1 

PVTZ/FS  CCSD+T*(3) 

FS  CCSD+T*(3) 

-6.9 

6.8 

a)  All  structures  were  optimized  using  the  DZP  basis. 
6.5  Relative  Energies 

Although  the  single  reference  and  Fock-space  methods  agree  on  the  geometries 
of  these  species,  they  do  not  agree  so  well  on  relative  energies.  As  shown  in  Ta- 
ble 6.2,  QRHF-CC  methods  and  FSCCSD  predict  the  2S+  structure  to  be  more  sta- 
ble by  varying  amounts,  while  FSCC  methods  including  triples  predict  the  2ng  state 
to  be  more  stable.  The  variations  in  the  FSCC  triples  results  due  to  choice  of  ac- 
tive space  are  on  the  order  of  0.2  kcal  mol"1  for  this  molecule,  so  they  cannot  have 
an  appreciable  effect  on  these  results.  There  are  two  important  trends  in  the  rela- 
tive energies,  one  related  to  the  FSCC  triples  methods,  and  the  other  due  to  the  ba- 
sis. 
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The  effect  of  increasing  the  basis  set  from  DZP  to  PVTZ  is  a  very  uniform  change 
of  approximately  2.5  kcal  mol1  in  the  relative  energy,  however  the  direction  of  the 
change  is  opposite  for  the  two  groups  of  methods:  For  the  QRHF-CC  and  FSCCSD 
methods,  the  gap  between  the  2ng  and  2S+  states  widens,  while  for  the  FSCCSD+T(3) 
and  FSCCSD+T*(3)  methods,  the  gap  is  reduced.  These  changes  match  what  would 
be  expected  if  the  change  in  basis  stabilized  the  2£+  state  while  leaving  the  2ng  state 
essentially  unchanged  for  all  methods.  Such  an  effect  could  result  if  the  2E+  state  has  a 
relatively  diffuse  character  which  for  which  the  change  to  a  PVTZ  basis  brings  significant 
improvement,  while  the  2ng  state  does  not  gain  any  special  advantage  from  the  larger 
basis. 

The  second  important  trend  in  is  the  marked  difference  in  the  relative  energies 
between  the  QRHF-CC  and  FSCCSD  methods  on  one  hand  and  FSCCSD+T(3)  and 
FSCCSD+T*(3)  on  the  other.  For  both  basis  sets,  there  is  change  of  10-15  kcal  mol"1 
in  the  relative  energies,  and  the  order  of  the  states  is  reversed. 

Looking  at  the  actual  energies  of  the  states,  rather  than  relative  energies,  reveals  that 
the  FSCC  third-order  triples  contribution  is  positive  for  the  £  states  and  negative  for 
the  n  state,  and  large  enough  to  change  the  order  of  the  states.  Comparison  with  the 
raw  data  from  the  previous  chapter  and  for  some  other  linear  carbon  clusters,  discussed 
below,  reveals  an  interesting  fact:  In  almost  every  case,  for  nonlinear  molecules,  the 
T(3)  and  T*(3)  corrections  to  FSCCSD  have  a  positive  sign,  regardless  of  the  symmetry 
of  the  state.  For  linear  molecules,  however,  the  pattern  is  that  £  states  have  a  positive 
triples  contribution,  while  II  states  have  a  negative  contribution.  This  means  that,  for  the 
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IP  calculations  on  linear  molecules,  triples  nearly  always  increase  the  IP  of  £  states  and 
decrease  that  of  II  states.  For  cases  in  which  the  relative  energy  is  of  interest,  the  results 
will  of  course  depend  on  the  relative  order  of  states  at  the  FSCCSD  level.  In  addition, 
nonlinear  structures  of  Cs+  calculated  in  the  course  of  this  work  (structures  with  no  more 
than  about  30°  deviation  from  linearity  in  any  one  angle),  the  same  pattern  is  obtained  for 
states  which  correlate  to  the  linear  S  or  n  states,  regardless  of  their  actual  symmetry.  For 
example,  in  Civ  symmetry,  the  first  two  states  of  Bi  symmetry  correlate  to  the  2ng  and 
2EU+  linear  states  respectively,  and  the  triple  excitation  effect  on  these  states  is  predicted 
by  which  state  they  correlate  to  in  the  linear  molecule. 

Such  a  pattern  is  not  inherent  in  the  FSCCSD+T(3)  or  FSCCSD+T*(3)  equa- 
tions, which  while  taking  into  account  the  symmetry  of  the  molecule,  in  no  way 
depend  on  it.  Computationally,  the  presence  of  counter-examples  among  the  vari- 
ous linear  molecules  surveyed  (one  II  state  of  02+  and  one  of  Cg"1")  confirms  that 
it  is  not  built  into  the  equations.  Also,  the  positive/negative  pattern  could  not  be 
established  for  nondegenerate/degenerate  states  molecules  in  degenerate  point  groups 
other  than  the  linear  ones.  It  would  seem  that  some  characteristic  of  the  elec- 
tronic structure  of  these  linear  molecules  causes  the  behavior,  probably  exagger- 
ated by  the  low  level  of  approximation  of  the  FSCC  triples  contribution.  The 
higher  order  terms  included  in  the  T*(3)  model  reduce  the  difference  in  ener- 
gies by  about  5  kcal  mol"1  compared  to  the  FSCCSD+T(3)  result,  so  improv- 
ing the  treatment  of  FSCC  triples  would  probably  produce  a  more  reasonable  an- 
swer. 


100 

Table  6.3.  Relative  energies  (in  kcal  mol"1)  of  B^^  C$+  stationary  points  for  several 
QRHF-CC  triples  approximations;  FSCC  values  are  also  given  for  reference. 

Relative  Energy 

Methoda)  .  . 
2ng  2EU+ 

QRHF-CCSD 
QRHF-CCSD[T] 
QRHF-CCSDT-lb 
QRHF-CCSDT-3 


-18.4 

0.0 

2.2 

0.0 

6.5 

0.0 

1.7 

0.0 

-6.1 

0.0 

-15.9 

0.0 

-13.0 

0.0 

FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 

a)  QRHF-CC  calculations  performed  at  the  QRHF-CCSD  optimized  geometry,  FSCC 
calculations  at  the  optimized  geometry  for  that  method.  All  calculations  used  a  DZP 
basis. 

The  question  that  must  now  be  answered  is  which  method(s)  provide  the  most  reliable 
relative  energies.  The  results  for  C3+,  for  both  QRHF-CC  and  FSCC  approaches,  as 
shown  in  Table  5.12,  do  not  provide  a  clear  answer,  but  some  useful  comparisons  can 
be  drawn.  In  the  single  reference  CC  case,  there  is  for  both  molecules,  a  significant 
change  in  the  relative  energy  upon  inclusion  of  the  perturbative  triples  correction.  This 
is  consistent  with  results  for  SRCC  method  applied  to  other  systems  that  would  normally 
be  considered  multireference  in  nature  —  the  CCSD  result  is  quite  misleading,  but  even 
the  lowest  order  addition  of  triple  excitations  gives  a  reasonable  result.  Going  to  the 
next  higher- level  inclusion  of  the  triple  excitation  effects,  CCSDT-1,  the  results  appear 
somewhat  worse  compared  with  other  triples  methods.  This  is  also  the  case  in  Cs+,  as 
shown  in  a  comparison  of  energies  at  the  two  Dod,  stationary  points  in  Table  6.3.     In 
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their  paper  on  C3+,  Watts  et  al.    [80]   attributed  this  behavior  to  the  large  fourth-order 
contribution  of  the  CCSD[T]  method  being  exaggerated  by  iteration  in  the  CCSDT-1 
method.  At  still  higher  levels,  the  triple  excitation  effect  stabilizes  for  both  C3+  and  Cs+. 

The  FSCC  results  have  also  been  shown  in  Table  6.3  with  the  same  zero  of  relative 
energy  as  for  the  QRHF-CC  data  in  order  to  simplify  comparison.  As  for  C3+,  the 
FSCCSD  result  is  not  in  line  with  the  QRHF-CC  numbers  including  triple  excitation 
effects.  For  Cs+,  the  difference  is  made  more  pronounced  by  the  fact  that  it  changes  the 
order  of  the  two  states.  The  inclusion  of  third-order  triples  effects  in  FSCCSD,  which 
has  already  been  shown  to  strongly  favor  the  n  state  over  the  S,  gives  relative  energies 
similar  to  the  QRHF-CCSD  single  reference  results.  The  same  pattern  is  seen  in  the  C3+ 
relative  energies  of  Table  5.12. 

Experience  with  QRHF-CC  methods  for  problems  of  this  sort  suggests  that,  when 
triples  excitation  effects  are  included,  the  results  are  quite  good,  and  in  this  case  we 
have  been  able  to  demonstrate  a  "convergence"  in  the  relative  energies  with  increasingly 
sophisticated  triples  approximations.  FSCC,  on  the  other  hand,  does  not  show  any 
convergence,  and  it  appears  that  the  E-TI  splitting  may  be  artificially  exaggerated  for 
linear  molecules  by  the  third-order  triples  approximations.  All  this  serves  to  make  the 
FSCC  relative  energies  less  believable  than  the  QRHF-CCSD[T]  ones.  By  combining 
the  trends  for  DZP/QRHF-CC  calculations  in  Table  6.3  with  the  data  in  Table  6.2,  the 
PVTZ/QRHF-CCSD[T]  relative  energies,  which  place  the  2ng  state  about  9  kcal  mol"1 
higher  than  the  2S+  state  should  be  good  predictor  of  results  from  higher-level  correlation 
treatments,  with  an  accuracy  of  ±4  kcal  mol"1.  Larger  basis  sets  would  be  expected  to 
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preferentially  stabilize  the  2S+  state,  because  of  its  more  diffuse  nature,  until  the  space 
is  saturated.  However  this  would  serve  to  increase  the  gap  between  the  states,  so  there 
is  little  likelihood  that  simply  enlarging  the  basis  would  produce  a  qualitative  change  in 
the  order  of  the  states. 

6.6  Stability  of  the  Stationary  Point  Structures 

To  this  point,  nothing  has  been  said  about  whether  or  not  these  structures  are  actually 
local  minima  or  some  other  type  of  stationary  point.  This  is,  of  course,  an  important 
question  because  if  these  structures  are  not  stable,  further  searching  will  be  required  to 
locate  true  local  minima.  This  is  especially  true  for  the  2ng  state,  which,  since  it  is 
degenerate,  is  susceptible  to  the  Renner-Teller  effect  (RTE). 

In  this  section,  harmonic  frequencies  are  presented  based  on  QRHF-CCSD,  FSCCSD, 
FSCCSD+T(3),  and  FSCCSD+T*(3)  calculations.  All  of  these  methods  have  problems 
in  these  circumstances,  which  will  be  discussed  as  appropriate.  However  it  should  be 
stated  up  front  that  QRHF-CCSD[T]  or  higher  level  single  reference  results  would  be 
most  desirable,  but  their  computational  cost  at  the  present  time  is  prohibitive  for  such  an 
application.  The  conclusions  drawn  here  have  been  checked  by  point-wise  exploration 
of  the  surfaces  away  from  linear  geometries. 

Table  6.4  shows  the  harmonic  frequencies  computed  for  the  2ng  structure.  All  four 
methods  agree  on  the  a  stretching  frequencies,  however  the  bending  modes  are  more 
complicated.  Both  QRHF-CCSD  and  FSCCSD  showed  some  lifting  of  the  degeneracy  of 
the  ix  modes,  as  might  be  expected  from  the  RTE.  The  FSCCSD+T(3)  and  FSCCSD+T*(3) 
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Table  6.4.  Harmonic  vibrational  frequencies  of  2ng  C$+  (in  cm"1)  using  various  theo- 
retical methods. 


Mode 

QRHF-CCSD 

FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 

IT 

137? 

187 

93i 

90i 

29 

43i 

■K 

88  i 

116 

21 

24 

55 

55 

TT 

242 

490 

494 

495 

496 

494 

a& 

760 

761 

757 

757 

<?u 

1262 

1306 

1292 

1288 

a& 

1828 

1895 

1885 

1880 

^u 

1884 

1904 

2121 

2023 

methods  did  not  produce  exactly  degenerate  modes,  but  the  differences  were  quite  small. 
The  primary  conclusion  from  these  data  and  from  additional  calculations  carried  out  at 
nonlinear  structures  is  that  the  PES  for  bending  from  this  state  is  extremely  flat.  The 
FSCC  data  required  using  numerical  differentiation  techniques  on  the  energy.  The  usual 
procedure  is  to  perturb  the  positions  of  the  nuclei  by  5  mBohr  to  obtain  the  necessary 
points.  In  this  case,  in  order  to  get  even  the  smallest  measurable  change  in  energy, 
the  step  size  for  the  bending  frequencies  had  to  be  increased  by  an  order  of  magnitude. 
But  with  such  a  large  step,  the  anharmonicity  of  the  surface  causes  problems.  As  a 
consequence,  the  FSCC  n  frequencies,  with  the  possible  exception  of  the  495  era"1 
mode,  should  not  be  given  any  credence.  The  QRHF-CCSD  frequencies  were  evaluated 
from  numerical  differentiation  of  analytic  gradients,  and  so  should  be  numerically  more 
reliable,  but  the  lack  of  triple  excitation  effects,  in  view  of  their  crucial  importance  for 
a  reasonable  representation  of  the  system,  make  the  QRHF-CCSD  bending  frequencies 
equally  questionable. 
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Table  6.5.  Harmonic  vibrational  frequencies  of  2£+  Cs+  (in  cm"1)  using  various  theo- 
retical methods. 


Mode 

QRHF-CCSD 

FSCCSD 

FSCCSD+T(3) 

FSCCSD+T*(3) 

7T 

148/ 

161 

nil 

94  i 

■K 

113 

273 

75 

94 

7T 

381 

478 

414 

411 

a 

771 

778 

790 

783 

a 

1452 

1453 

1498 

1474 

a 

2100 

2142 

2054 

2093 

a 

2237 

2291 

2238 

2261 

For  2E+  C5+,  the  vibrational  frequencies  shown  in  Table  6.5  are  somewhat  better 
than  those  for  the  2ng  state.  Once  again,  the  stretching  frequencies  are  in  accord  for  all 
four  methods.  For  the  n  modes,  the  surface  is  fairly  flat  compared  to  most  molecules 
(but  less  so  than  for  the  2ng  state),  so  numerically,  the  FSCC  results  have  a  sizeable 
margin  of  error.  However  for  the  same  reasons  cited  above,  the  QRHF-CCSD  results 
should  not  be  entirely  trusted,  leaving  FSCCSD  as  the  best  choice.  As  a  check  for  the 
nature  of  a  possible  nonlinear  minimum,  the  FSCCSD+T*(3)  minimum  was  located.  This 
is  a  A  state,  which  has  one  terminal  carbon  bent  away  Irom  the  linear  axis  by  about 
10°.  At  the  DZP/FSCCSD+T*(3)  level,  this  geometry  is  about  4  kcal  mol"1  more  stable 
than  the  linear  2E+  state.  The  small  energy  difference  attests  to  the  anharmonicity  of 
the  surface  in  this  region. 

Thus  vibrational  frequency  calculations  indicate  that  both  states  are  probably  linear, 
but  floppy  —  particularly  the  2ng  state.  Even  if,  with  higher-level  correlated  methods, 
it  is  found  that  there  are  in  fact  nonlinear  minima,  there  is  no  evidence  to  suggest  that 
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the  ordering  of  the  states  will  change.  Due  to  the  very  large  anharmonicities  of  bending 
modes  on  both  surfaces,  the  computed  n  frequencies  should  not  be  considered  reliable. 
However  the  stretching  frequencies  should  be  good,  and  may  be  useful  in  experimental 
studies  of  this  molecule. 

6.7  Implications  for  Other  Carbon  Clusters 

The  numerous  complexities  encountered  for  Cs+  prompt  questions  about  what  the 
situation  might  be  for  other  odd-numbered  cationic  carbon  clusters.  In  a  sense,  the 
"problems"  of  Cs+  stem  from  the  fact  that  it  has  so  many  states  close  in  both  energy  and 
geometry,  in  such  a  way  that  they  can  interact  strongly.  Using  the  Fock-space  coupled 
cluster  approach,  it  is  relatively  easy  to  explore  the  structure  of  other  clusters  and  see 
how  their  electronic  structure  compared  to  that  of  Cs+. 

For  these  calculations,  the  6-31G*/SCF  optimized  geometries  of  the  neutral  clusters 
of  Raghavachari  and  Binkley  [96]  were  used,  since  they  provide  a  consistent  set  of 
structures  though  C9.  DZP/FSCCSD  calculations  were  performed  at  these  geometries  to 
obtain  the  vertical  energy  differences  between  the  low-lying  states.  The  results  of  these 
calculations  are  shown  in  Table  6.6. 

There  are  two  importance  effects  to  note  in  this  data.  First,  the  relative  energy  of  the 
two  £  states  (highlighted  in  the  table)  decreases  from  C3+  to  a  crossover  point  at  Cj+, 
where  they  are  essentially  degenerate,  then  in  C9+  the  gap  begins  to  widen  again,  but 
with  the  gerade  symmetry  state  lower.  The  second  effect  is  that  of  the  lowest  II  state, 
which  alternates  between  gerade  and  ungerade,  and  drops  from  18  kcal  mol*1  above  the 
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Table  6.6.  Relative  energies  of  low-lying  electronic  states  of  carbon  cations  (in  kcal 
mol"1). 

Molecule  State  and  relative  energy 


C3+ 


C5+ 


c7+ 


c9+ 


2T  +  2Y  ■  2TT 

£»u  Zjg  llg 

0  18.77  36.44 


2ju  2jg  Hu  llg 

0  1.71  2.67  69.88 


2ng  2su+  2sg-  2nu 

0  11.91  11.93  57.29 


2n  2v  ■  2v  +  2TT 

llu  .Ug  .iju  llg 

0  20.83  20.95  50.47 


E  states  to  21  kcal  mol"1  below  them.  These  two  effects  combine  in  Cs+  with  that  the 
result  that  three  states  are  (at  this  geometry)  all  within  3  kcal  mol"1. 

Based  on  these  data,  the  pseudo-Jahn-Teller  interaction  of  the  two  S  states  would 
peak  for  C7+,  where  they  are  essentially  degenerate,  then  decrease  as  the  cluster  lengthens 
and  the  states  separate  again.  As  the  cluster  lengthens,  we  would  also  expect  the  energy 
lowering  from  the  D^  symmetry  to  the  CooV  minimum  to  decrease,  since  small  changes 
in  bond  length  would  have  a  lesser  effect  on  the  energy  of  the  overall  system.  At  the 
same  time,  gap  between  the  n  state  and  the  Ss  widens,  making  it  unlikely  that  even  the 
distorted  CooV  £  state  would  be  competitive  for  the  ground  state.  Finally,  the  2U  state 
is  subject  to  the  Renner-Teller  effect,  which  would  likely  (but  not  necessarily)  result  in 
a  nonlinear  ground  state. 
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Calculations  at  the  optimized  2SU+  structure  of  C3+,  given  in  the  previous  chapter, 
and  at  various  Cs+  stationary  points  shown  above  indicate  that  the  somewhat  arbitrary 
choice  of  geometry  used  in  the  calculations  presented  in  Table  6.6  does  not  affect  the 
picture  it  portrays,  merely  changing  the  exact  relative  energies  somewhat. 

Experimentally,  two  distinct  isomers  of  both  Ci+  and  C94"  have  been  observed  [84,88] 
with  differing  mobilities  or  chemical  reactivities.  These  have  been  referred  to  as  "linear" 
(reactive)  and  "cyclic"  (unreactive),  but  the  basis  for  these  structural  assignments  is  not 
strong.  Parent  and  McElvany  [84]  relied  on  very  simple  calculations  such  as  those  of 
Bernholc  and  Phillips  [97]  at  the  semiempirical  MNDO  level  which  considered  a  very 
restricted  set  of  possible  structures.  Von  Helden  and  coworkers  [88]  also  measured  the 
mobilities  of  linear  and  cyclic  alkanes  and  found  that  their  ratios  matched  those  observed 
for  the  two  isomers  of  the  pure  carbon  species.  In  fact,  since  "linear"  alkanes  are  not 
actually  linear  in  the  strict  sense  usually  applied  to  pure  carbon  clusters,  von  Helden's 
experiments  suggest  that  the  difference  between  a  truly  linear  cluster,  and  one  that  is 
distorted  from  linearity  could  not  be  resolved.  Although  this  may  be  viewed  as  a  matter 
of  semantics,  application  of  the  term  "linear"  in  this  context  is  somewhat  misleading; 
it  is  likely  that  the  reactive  "linear"  isomer  seen  in  experiments  is,  in  fact,  nonlinear 
due  to  the  Renner-Teller  effect.  Spectroscopic  or  Coulomb  explosion  studies  on  these 
cations  will  be  required,  probably  in  conjunction  with  more  detailed  theoretical  studies, 
to  resolve  this  question. 
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6.8  Recapitulation 

The  ground  state  of  Cs+  has  been  found  to  be  2E+,  with  Coov  symmetry,  with  a  2ng 
state  9  kcal  mol"1  higher  in  energy.  Harmonic  stretching  frequencies  of  the  two  structures 
have  been  computed  and  should  help  spectroscopists  in  identifying  this  molecule.  With 
respect  to  bending  motions,  both  surfaces  are  rather  flat  —  especially  the  2ng  state.  At 
the  FSCCSD  level,  the  2£+  state  is  stable,  but  with  other  methods  it  is  a  transition  state. 
However  even  if  this  state  is  found  to  be  nonlinear  with  more  reliable  methods,  the  change 
in  energy  would  be  expected  to  be  small,  and  would  not  change  the  important  conclusions. 

Although  they  were  not  deemed  the  most  reliable  results  forjudging  relative  energies, 
etc.,  FSCC  methods  proved  invaluable  in  this  study.  They  provided  the  first  insight  into 
the  nature  of  the  problem  that  had  plagued  other  calculations,  and  were  used  to  explore 
several  PES  at  once  for  interesting  features.  The  third-order  triples  approximations, 
FSCCSD+T(3)  and  FSCCSD+T*(3),  appear  suffer  from  an  interaction  between  the  level 
of  approximation  and  the  electronic  structure  of  linear  molecules  that  uniformly  favors 
n  states  by  a  substantial  amount.  It  is  expected  that  better  treatment  of  the  triples 
contribution  will  produce  more  reasonable  results. 

FSCC  methods  have  also  been  applied  to  other  clusters  in  the  Cn+  (n  odd)  series, 
and  it  appears  that  Cs+  is  the  worst  possible  case  with  respect  to  proximity  of  the  states. 
As  the  size  of  the  cluster  increases,  the  2n  state  drops  in  energy  with  respect  to  the  2E, 
so  that  for  largely  molecules,  the  ground  state  almost  certainly  derives  from  the  2IL  Of 
course  since  it  is  a  degenerate  state,  Renner-Teller  effects  may  result  in  the  actual  ground 
state  nonlinear,  but  this  will  require  a  significant  computational  effort  to  establish. 


CHAPTER  7 
SUMMARY 


The  Fock-space  coupled  cluster  theory  has  been  derived  and  its  relationship  with 
the  more  familiar  single  reference  CC  theory  discussed.  The  full  FSCCSDT  equa- 
tions for  the  (0,1)  sector  (one  electron  ionizations)  have  been  presented  in  the  form 
of  spin  orbital  equations  along  with  the  necessary  elements  of  the  coupled  cluster  ef- 
fective Hamiltonian,  //jv,  also  in  spin  orbital  form.  Several  perturbative  and  itera- 
tive approximations  to  the  full  FSCCSDT  model  are  proposed  as  a  compromise  be- 
tween the  cost  of  the  complete  inclusion  of  triple  excitation  effects  and  the  improved 
accuracy  they  should  afford.  The  approximations,  as  presented,  do  not  have  certain 
invariance  properties  which  might  be  desirable.  In  SRCC  methods,  experience  has 
shown  that  the  incorporation  of  the  terms  required  to  obtain  the  desired  invariance 
generally  leads  to  improved  performance  of  the  approximate  triple  excitation  method. 
The  modifications  to  the  proposed  FSCC  triples  approximations  to  meet  these  require- 
ments were  discussed.  It  was  also  observed  that  approximating  the  triple  excitation 
contribution  causes  the  method  to  lose  the  invariance  with  respect  the  choice  of  ac- 
tive space  which  would  apply  to  either  the  FSCCSD  or  FSCCSDT  methods.  The 
proof  of  this  invariance  relies  on  satisfaction  of  the  Fock-space  Bloch  equation,  so  it 
is  not  clear  that  any  approximate  method  can  be  formulated  with  the  property,  and  the 
severity  of  the  problem  in  practice  will  have  to  be  determined  by  computational  exper- 
iments. 
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As  a  first  step  in  understanding  triple  excitation  effects  in  FSCC  methods,  two  of  the 
simplest  approximations  have  been  implemented  for  use  with  both  closed-  and  open-shell 
reference  functions.  The  FSCCSD+T(3)  method  includes  only  third-order  contributions 
of  triples  to  the  energy,  and  the  FSCCSD+T*(3)  method  incorporates  some  higher-order 
effects  by  replacing  raw  H^  elements  used  in  the  evaluation  of  the  T(3)  contribution 
with  Hn.  The  spin  integrated  equations  from  which  the  computer  program  was  written 
are  given  in  the  Appendix.  These  methods  were  then  applied  to  a  variety  of  molecules  in 
order  to  get  a  general  feeling  for  how  the  FSCC  third-order  triples  approximation  perform 
compared  to  FSCCSD  and  single  reference  methods. 

Computational  experiments  with  different  active  spaces  show  that  in  most  cases,  the 
uncertainty  introduced  is  small.  There  are,  however,  a  few  molecules  for  which  the  effect 
is  more  pronounced.  Ideally,  it  would  be  useful  to  be  able  to  predict  which  systems  will 
have  the  largest  uncertainty,  but  no  common  thread  could  be  found  among  the  cases 
studied;  more  experience  with  FSCC  triples  approximations  will  be  required  to  reach  that 
level.  It  is,  however,  a  fairly  simple  matter  to  test  the  effect  on  any  given  molecule  of 
interest.  Calculations  of  ionization  potentials,  mostly  using  DZP-quality  basis  sets  show 
that  FSCCSD  and  FSCCSD+T*(3)  both  perform  about  as  well  as  QRHF-CCSD[T]  in 
comparison  with  experiment,  with  an  average  absolute  error  of  about  0.2  eV.  The  higher- 
order  effects  included  in  the  T*(3)  model  appear  to  have  a  moderating  effect  on  the  pure 
third-order  correction,  and  thus  compare  better  with  experimental  results.  Fock-space 
methods  should  also  be  well-suited  to  the  study  of  systems  exhibiting  symmetry  breaking 
or  other  vibronic  interaction  effects,  and  several  such  systems  are  examined.  Although 
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there  is  no  conclusive  experimental  or  theoretical  data  to  compare  with,  it  appears  that 
FSCC  methods  are  quite  capable  of  obtaining  good  geometries,  and  in  most  cases  good 
vibrational  frequencies.  For  relative  energies,  however,  the  results  are  not  so  satisfying, 
either  with  the  basic  FSCCSD  method,  or  when  it  is  augmented  with  triples  —  while 
SRCC  methods  tend  to  show  a  "convergence"  in  the  relative  energy  with  increasingly 
sophisticated  treatment  of  triples  excitations,  FSCCSD  relative  energies  begin  several 
kcal  mol"1  away  from  the  SRCC  results  and  are  not  obviously  improved  by  third-order 
triples.  In  fact  throughout  these  calculations,  there  were  signs  of  the  basic  differences 
between  the  SRCC  excitation  operator  hierarchy  and  the  FSCC  hierarchy,  which  were 
described  from  a  theoretical  standpoint  earlier. 

A  combination  of  SRCC  and  FSCC  methods  was  applied  to  the  problem  of  lo- 
cating the  ground  state  structure  of  Cs+,  attempting  to  use  the  advantages  of  each 
method.  The  FSCC  methods  proved  invaluable  for  "exploring"  the  potential  energy 
surface  and  focusing  on  regions  of  interest  or  uncovering  potential  problems  such  as 
strong  vibronic  interactions.  In  many  electronic  structure  problems,  regions  of  degen- 
eracy or  near-degeneracy  are  known  a  priori  not  to  provide  the  sought-after  property 
(i.e.  local  minimum),  and  by  recognizing  them,  computational  headaches  can  often  be 
avoided. 

During  the  course  of  the  Cs+  calculations,  it  became  apparent  that  the  FSCCSD+T(3) 
and  FSCCSD+T*(3)  methods  show  a  very  strong  pattern  of  triples  positive  correc- 
tions for  £  states  and  negative  for  n  states,  thus  energetically  favoring  the  n  state. 
This  pattern  is  not  observed  in  general,  nonlinear  molecules.     It  is  concluded  that 
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these  effects  are  probably  due  to  some  aspect  of  the  electronic  structure  of  linear 
molecules  being  exaggerated  by  the  low-order  approximation  used  for  triple  excita- 
tions. 

Despite  these  problems,  the  ground  state  of  Cs+  was  determined  to  be  2  E+, 
with  Coov  symmetry,  with  a  2ng  state  roughly  9  kcal  mol"1  higher  in  energy.  Both 
of  these  states  appear  to  be  linear,  but  have  relatively  flat  potential  surfaces  with 
respect  to  bending.  Fock-space  coupled  cluster  calculations  have  also  been  performed 
on  the  series  of  molecules  C„+,  (n=3,5,7,9)  in  order  to  examine  the  implications  of 
the  C5+  results  for  the  rest  of  the  series.  While  problems  with  pseudo-Jahn-Teller 
interaction  between  the  2Sg"  and  2EU+  states  grows  worse  as  the  cluster  gets  larger, 
the  n  state  becomes  much  more  stable  than  the  Es  in  Cn+  and  Cg*.  Thus  C5+  is 
perhaps  the  worst  of  these  clusters  with  respect  to  the  complexity  of  its  electronic 
structure,  and  the  lowest  linear  energy  linear  structures  of  C-j+  and  C9+  are  predicted 
to  be  2n,  with  the  possibility  of  distortion  from  strict  linearity  due  to  the  Renner-Teller 
effect. 

There  is  one  further  aspect  of  these  calculations  which  was  mentioned  very  little 
in  the  previous  chapters  —  the  cost  of  the  FSCC  calculations.  In  the  most  common 
application  of  FSCC  (0,1)  methods,  the  (0,0)  sector  reference  is  a  well-behaved  closed- 
shell  molecule.  Consequently,  the  entire  (0,0)  calculation  can  be  done  using  an  RHF 
reference.  Single  reference  CC,  on  the  other  hand,  would  require  an  open-shell  reference 
function,  and  therefore  an  open-shell  CC  calculation,  which  is  about  four  times  more 
expensive  than  for  a  closed-shell  reference  (the  FSCC  calculation  itself  is  a  small  fraction 


113 
of  the  cost  of  the  reference  CC  calculation).  In  addition,  the  SRCC  calculation  is  much 
more  susceptible  to  convergence  problems,  so  the  calculations  may  take  more  iterations  of 
the  CC  equations  than  the  FSCC  (0,0)  closed-shell  calculation.  Although  FSCC  methods 
applied  to  open-shell  reference  functions  don't  have  this  advantage  over  SRCC,  they 
do  retain  the  second  advantage  of  the  Fock-space  approach  —  the  ability  to  evaluate 
several  states  in  a  single  calculation.  Single  reference  CC  calculations  of  vertical  IPs  or 
excitation  energies  require  one  calculation  per  state  involved,  while  FSCC  can  provide 
them  all  in  a  single  calculation. 

The  computational  advantages  of  the  FSCC  approach  mean  that  a  much  more  thor- 
ough examination  of  some  complex  electronic  structure  problems  is  possible  compared  to 
exclusive  use  of  single  reference  methods,  and  although  currently  available  FSCC  meth- 
ods do  not  seem  to  perform  as  well  as  might  be  desired  for  some  aspects  of  the  problem, 
this  will  often  be  outweighed  by  being  able  to  explore  the  system  better.  As  a  result,  it  is 
well  worth  putting  additional  effort  into  improving  our  understanding  of  FSCC  methods 
and  their  computational  performance.  With  respect  to  the  triples  approximations  investi- 
gated, it  is  clear  that  the  excitation  operator  hierarchy  of  SRCC  and  FSCC  is  different, 
both  theoretically  and  computationally.  While  low-order  triples  approximations  perform 
quite  well  in  most  cases  for  SRCC  methods,  they  are  not  so  successful  for  FSCC.  The 
obvious  next  step  is  to  move  to  higher-level  approximations,  and  look  for  better  perfor- 
mance. But  from  the  development  of  single  reference  triples  approximations,  it  would 
seem  that  the  most  useful  tool  is  a  large  base  of  calculations  which  make  it  possible  to 
identify  characteristics  of  molecules  that  lead  to  particularly  good  or  poor  results  from 
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approximate  triples  methods,  thereby  better  understanding  existing  methods  and  allow- 
ing more  refined  approach  to  formulating  better  approximations.   Hopefully,  this  work 
represents  the  beginning  of  such  an  effort  for  the  Fock-space  coupled  cluster  method. 


APPENDIX 
IMPLEMENTATION  OF  THE  PERTURBATIVE  TRIPLES  MODELS 


A.l   Introduction 

The  various  approximate  triple  excitation  models  presented  in  this  work  have  been 
implemented  in  the  ACES  II  program  system  [98]  by  building  upon  the  FSCCSD  program 
written  by  John  F.  Stanton.  Due  to  intrinsic  differences  of  these  methods  from  the 
majority  of  other  which  have  been  implemented  within  the  ACES  II  framework,  and 
also  due  to  differences  among  the  programmers,  there  are  some  aspects  of  the  FS  triples 
implementation  which  bear  comment  in  order  to  assist  anyone  who  may  wish  to  build  on 
this  work  in  the  future.  In  addition,  the  FSCC  literature  in  particular  tends  to  avoid  giving 
explicit  equations,  preferring  the  more  compact  operator  or  diagrammatic  expressions. 
This  makes  for  less  cumbersome  papers,  but  can  lead  to  ambiguity  as  to  exactly  which 
terms  are  present  in  a  given  method.  Consequently,  I  will  present  both  the  spin  orbital 
equations  and  the  final  spin-integrated  (suitable  for  general  reference  functions)  which 
have  been  coded. 

A. 2  Spin  Orbital  Equations  for  the  T(3)  and  T*(3)  Models 

The  second-order  triple  excitation  amplitudes  of  the  T(3)  model  consist  of  four  terms, 
two  from  III n,  the  three-body  part  of  Hn,  and  two  from  the  contraction  of  converged 
Tlj  '     with  Wtf.  The  restriction  on  order  means  that  only  the  Wff  component  of  Wn 
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will  actually  contribute. 

DiaHj  =  ~  (-)pP(ki/J)(im\\ki)t%  (A.la) 

+  (-yP(k/ij)(la\\keytbj  (A.lb) 

-(-)pP(k/ij,ab)(mb\\ij)rlam  (A.lc) 

+  (-)pP(ki/j)(ab\\ej)Tti.  (A.ld) 

The  notation  used  above  is  that  amplitudes  from  T'0'0'  are  denoted  with  t,  and  those 
from  J1'0,1)  with  r;  the  dots,  /  are  used  to  indicate  that  /  is  restricted  to  be  an  active 
orbital,  by  analogy  with  the  double  arrow  used  to  indicate  active  lines  in  diagrams.  D£~ 
is  the  coupled  cluster  denominator, 


D[fj   =  fkk  +  fii  +  fn  -  ff,  -  faa  -  fbb-  (A.2) 


Remember  that  throughout  this  appendix,  the  target  amplitude  rjf-j  is  only  correct  through 
second  order.  Additional  notation  to  indicate  this  would  be  too  cumbersome. 

The  triples  amplitudes  contribute  to  the  energy  via  the  effective  Hamiltonian, 

{HN,e«)rk  =  \(ij\\ab)ij#.  (A3) 


To  obtain  the  T*(3)  model  proposed  by  Pal,  Rittby,  and  Bartlett  [56],  two  modifica- 
tions to  the  above  equations  are  required.  First,  all  H^  integrals  appearing  in  Eqs.  A.l, 
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A. 2,  and  A. 3  are  replaced  by  the  corresponding  H^  elements: 

Jpq  >  Jpq 

(A.4) 

{pq\\rs)  ►  {pq\\rs) 

Second,  the  denominator  is  modified  to  incorporate  some  effects  from  the  renormalization 

term,  J1]0'1^  <—  H^^T^  '  \  by  replacing  the  /  term  with  the  corresponding  diagonal  of 

H^e{{.   The  resulting  equations  are: 

D^f  =  -(-fP(ki/j)(lm\\ki)t^ 

+  (-)pP(k/rj)(ia\\ke)tf 

(A.5) 

-(-YP(k/ij,ab)(mb\\ij)rlam 

+  (-)PP(ki/j)(ab\\ej)rtl, 

Dif  =  fkk  +  fa  +  fa  -  {HN>eK)n  -  faa  -  fa,  (A.6) 


and 


(HN^)ii  =  \{ij\\ab)r^.  (A.7) 


It  is  worth  observing,  however,  that  in  H^^,  the  integral  {ij\\ab)  and  the  Hn  element 
{ij\\ab)  are  identical,  so  that  Eqs.  A.l  and  A.5  are  in  fact  the  same  except  for  the  TJ  ' 
amplitude. 

A.3  Details  of  the  Implementation 

Implementation  of  the  T(3)  and  T*(3)  model  in  the  ACES  II  program  system 
was  done  both  with  an  eye  to  the  special  features  of  the  third-order  models,  and 
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at  the  same  time  providing  a  base  for  higher- level  perturbative  and  iterative  mod- 
els up  to  full  FSCCSDT.  Also,  as  this  is  the  first  implementation  of  these  methods 
within  the  ACES  II  system,  performance  concerns  were  subordinated  to  clarity  and 
generality  of  the  code.  As  a  result,  it  may  be  possible  to  improve  the  performance 
of  the  third-order  methods  by  specializing  the  code  to  them,  however  I  think  it  is 
worth  waiting  on  this  until  there  is  a  larger  body  of  experience  with  FS  triples  meth- 
ods. 

In  order  to  avoid  storing  the  complete  set  of  T3  '  amplitudes,  their  calculation 
is  organized  as  a  loop  over  the  required  (k,l)  pairs,  and  at  each  iteration  evaluating 
contribution  of  a  subset  of  TJ  '  '  having  the  same  size  as  T.?  '  to  i//v,eff-  This  approach 
should  work  for  all  methods  not  requiring  a  T3l  <—  T\  '  contribution.  Fixing  (  kj\ 
also  requires  the  ability  to  read  integral  and  amplitude  lists  with  these  indices  fixed.  The 
original  FSCCSD  implementation  includes  a  set  of  routines  which  has  come  to  be  known 
as  the  "Fock  space  1/0"  (FSIO)  package,  which  allows  the  use  of  subsets  with  indices 
restricted  to  active  or  inactive  via  scatter  and  gather  operations  using  precomputed  index 
vectors  for  each  subset.  This  package  was  extended  to  allow  for  "temporary  I/O"  (TIO) 
vectors  which  can  be  easily  changed  to  correspond  to  each  (k,l  J  in  the  loop.  This 
avoids  having  to  store  precomputed  scatter/gather  vectors  for  all  of  the  (k,l)  at  a  very 
small  cost  in  CPU  time.  For  methods  that  require  storage  of  the  full  TJ  '  '  from  one 
iteration  to  the  next,  some  changes  will  have  to  be  made.  Since  the  current  code  is 
structured  around  the  ability  to  store  objects  the  size  of  T]  '  in  memory,  it  should  be 
fairly  easy  to  store  the  T\  '     amplitudes  so  that  each  record  ("distribution"  in  ACES  II 
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terminology)  holds  an  entire  T2l  '  -sized  quantity  corresponding  to  a  given  (k,l),  in 
the  same  fashion  as  is  presently  done  in  the  SRCC  triples  codes  which  require  storage 
of  the  amplitudes. 

Probably  the  most  distinctive  aspect  of  this  implementation  is  the  treatment  of  the 
spin  cases  of  T3  in  deriving  the  final  form  of  the  equations.  Normally,  using  the 
permutational  symmetry  of  the  amplitudes,  one  would  consider  four  unique  spin  cases. 
In  order  to  simplify  evaluation  of  approximations  which  don't  need  the  full  //^,eff  matrix 
(spectator  and  active-active  subsets,  for  example),  the  indices  ( k,  l)  are  treated  as  distinct 
from  the  remaining  indices,  (i,j,a,b).  This  leads  to  six  distinct  spin  cases  —  two  from 
the  (k,l)  part  times  three  from  the  other  four  indices  —  but  keeps  the  code  very  simple 
with  respect  to  the  use  of  RHF  vs.  UHF  references  and  the  contribution  of  the  amplitudes 
to  Hfteft.  If  the  complete  r3  '  tensor  must  be  evaluated,  as  for  higher-level  methods,  it 
would  make  sense  to  treat  the  minimal  set  of  spin  cases  and  make  full  use  of  permutational 
symmetry.  This  will  primarily  involve  reorganization  of  the  contraction  of  the  amplitudes 
into  //jvieff,  and  of  course  avoiding  of  the  redundant  terms  in  generating  the  amplitudes. 

Because  of  the  obvious  similarity  of  the  T(3)  and  T*(3)  models,  they  have  been 
implemented  using  a  single  set  of  routines  capable  of  either  method  rather  than  two  sets 
of  routines  distinguished  only  by  accessing  slightly  different  integral  lists.  This  is  done 
by  writing  the  code  for  the  T*(3)  method  and  using  an  array  to  change  the  list  numbers 
from  those  for  H^  to  those  for  H^.  (T*(3)  is  chosen  as  the  "lead"  method  because  the 
Hn  integrals  it  requires  lack  the  bra-ket  transposition  symmetry  of  the  H^  integrals,  so  a 
code  that  is  correct  for  T*(3)  will  be  correct  for  T(3)  when  the  integrals  are  changed,  but 
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generally  not  for  the  reverse.)  The  same  mapping  of  list  numbers  is  used  to  handle  the  fact 
the  for  closed  shell  reference  functions,  many  lists  become  redundant  and  are  not  stored. 
In  most  other  ACES  II  code,  different  cases  are  used  within  the  code  for  UHF  vs.  RHF 
references;  in  the  FS  triples,  the  approach  is  to  use  the  same  code  with  the  list  numbers 
properly  mapped,  and  limit,  at  a  higher  level  in  the  code,  which  spin  cases  are  evaluated. 
It  must  also  be  noted  that  the  evaluation  of  H^  in  the  program  LAMBDA  overwrites 
some  of  the  pure  H^  integrals,  so  it  is  necessary  to  preserve  them  in  order  to  run  an 
FSCCSD+T(3)  calculation  later.  This  is  done  with  a  very  short  program,  SAVEINTS, 
which  must  be  run  before  LAMBDA,  and  simply  copies  the  relevant  integrals  to  other 
(unused)  list  numbers.  The  list  mapping  mechanism  then  simply  redirects  references  to 
the  integrals  lists  to  their  new  location  rather  than  copying  the  integrals  back  into  their 
original  locations.  By  mapping  the  lists  instead  of  copying  them,  it  is  possible  to  do  both 
T(3)  and  T*(3)  calculations  after  a  FSCCSD  in  a  single  invocation  of  the  program. 

A.4  Final  Equations  for  the  T(3)  and  T*(3)  Models 

In  Tables  A.l — A.4  the  final  spin  integrated  form  of  the  equations  used  in  the  third- 
order  methods,  suitable  for  a  general  open  shell  reference,  is  displayed.  Six  spin  cases  are 
shown,  and  the  resulting  terms  are  grouped  according  to  the  natural  pairing  of  indices  for 
the  loop  over  irreducible  representations  of  the  TJ  '  — like  part  of  the  amplitude.  Only 
the  T(3)  equations  are  shown,  but  as  discussed  above,  the  T*(3)  equations  are  identical 
except  that  H^  integrals  are  substituted  for  H^. 
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Table  A.l.  Final  spin  integrated  equations  for  Eq.  A. la  (also  shown  at  the  top  of  the 
table).  Uppercase  indices  denote  alpha  spin  orbitals,  lowercase  beta  spin.  The  symbol  P 
in  the  equations  denotes  the  permutation  operator  (—)pP(ij). 


-(-yp(ki 

fm 

Hl«)*ft 

Spin  Case 

Loop  over  T^ 

AB  IJ  KL 

+P{LM\\KI)tjB 

+  (LM\\U)tAKBM 

ab  ij  KL 

+  P(Lm\Ki)tabn 

Ab  Ij  KL 

-{LM\\Kl)tH 

-(Lm\Kj)tfbn 

+(Lm\Ij)tfm 

AB  IJ  kl 

+P{Ml\Ik)tjB 

ab  ij  kl 

+P(lm\\ki)tfm 

+(lrn\\ij)tll 

Ab  Ij  kl 

-(Ml\Ik)t^\ 

-{lm\\kj)tfl 

+(Ml\Ij)t$k 

Table  A. 2.  Final  spin  integrated  equations  for  Eq.  A. lb  (also  shown  at  the  top  of  the 
table).  Uppercase  indices  denote  alpha  spin  orbitals,  lowercase  beta  spin.  The  symbol  P 
in  the  equations  denotes  the  permutation  operator  {  —  )pP(ab). 


-(-yp(k/ij, 

ab)(la\\ke)t\l 

Spin  Case 

Loop  over  TtJ 

Loop  over  I\a 

Loop  over  r,-j 

AB  IJ  KL 

-P(LA\\KE)tff 

+P{LA\\IE)tj» 

-P(LA\\JE)tfg 

ab  ij  KL 

-P(La\Ke)th* 

+P(La\Ei)tf1 

-P(La\Ej)tft 

Ab  Ij  KL 

+(LA\\KE)tfb 
+(Lb\Ke)tf? 

-(LAUIE^fj 
-(Lb\Ej)tfg 

-(Lb\Ie)4) 

AB  IJ  kl 

-P{lA\kE)tff 

+  P(lA\eI)tfk 

-P(lA\eJ)t% 

ab  ij  kl 

-P(la\\ke)t% 

+P(la\\ie)lfi 

-P(la\\je)tfk 

Ab  Ij  kl 

+  {lA\kE)tfjb 
+(lb\\ke)tff 

-(lA\eI)tlb3 
-(lb\\je)tf£ 

-mjE)tfkb 

A. 5  Testing 


Of  course  programs  are  rarely  exactly  correct  when  they  are  first  written,  and  they 
must  be  tested  to  verify  their  correctness  and  point  up  bugs.  This  section  lists  some  of 
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Table  A. 3.  Final  spin  integrated  equations  for  Eq.  A.lc  (also  shown  at  the  top  of  the 
table).  Uppercase  indices  denote  alpha  spin  orbitals,  lowercase  beta  spin.  The  symbol  P 
in  the  equations  denotes  the  permutation  operator  (  —  )pP(ab). 


-(-)pP(k/ij,ab)(mb\\ij)T 


la 


Spin  Case 


Loop  over  T 


i  j 


Loop  over  Tt 


Loop  over  T^ 


AB  IJ  KL 


LA 


.LA 


P(MB\\IJ)tJ&       -P(MB\\KJ)t^j      +P(MB\\KI)t^j 


ab  ij  KL 


-P(mb\\ij)TJ! 


La 
m 


La 


-P(Mb\Kj)Tti 


+P{Mb\Ki)rLMa1 


Ab  Ij  KL 


iMb\Ij)r^AM 
(Am\Ij)4 


LA 


■(Mb\KJ)rjfi 


Lb 
m 


Lb 


-{MA\\KI)rti 


+(Am\Kj)rj 


Lb 
m 


AB  IJ  kl 


-P(MB\\IJ)Tfik         -P{Bm\Jk)r^  +  P(Bm\Ik)rj 


Al 

in 


.Al 
Jm 


ab  ii  kl 


■P(mb\\ij)T'k 


la 


■P{mb\\kj) 


.la 


+  P(mb\\ki)T^ 


Ab  Ij  kl 


■(Mb\Ij)T*lk 
■{Am\Ij): 


km 


■{mb\\kj)ri 

(Am\Ik) 


Al 
m 

.lb 
mi 


+  {Mb\Ik)T^ 


Table  A.4.  Final  spin  integrated  equations  for  Eq.  A. Id  (also  shown  at  the  top  of  the 
table).  Uppercase  indices  denote  alpha  spin  orbitals,  lowercase  beta  spin.  The  symbol  P 
in  the  equations  denotes  the  permutation  operator  (—)pP(ij). 


(-)pP(k, 

/.?)<< 

'6l|ej>/.- 

Spin  Case 

Loop  over  T^ 

AB  IJ  KL 

+P(AB\\EJ)4f 

+(AB\\EK)iff 

ab  ij  KL 

+P{ab\\ej)TJfi 

Ab  Ij  KL 

+  (Ab\Ej)rtf 

+{Ab\Ie)r% 

-{Ab\Ke)rf* 

AB  IJ  kl 

+P(AB\\EJ)tF< 

ab  ij  kl 

+P(ab\\ej)rj;i 

+{a.b\\ek)T% 

Ab  Ij  kl 

+  (Ab\Ej)T<f 

+{Ab\Ie)rti 

-[Ab\Ek)rfj 

the  tests  that  have  been  applied  to  the  FS  triples  code  in  the  course  of  its  development. 
While  they  are  probably  obvious  to  most  people,  all  have  proven  useful  and  none  should 
be  forgotten. 
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1.  Symmetry.  Given  a  molecule  with  point  group  symmetry,  the  code  should  give  the 
same  result  when  run  in  C,  or  in  full  symmetry. 

2.  UHF  vs.  RHF.  For  a  closed  shell  molecule,  using  UHF  and  RHF  reference  functions 
should  give  the  same  result. 

3.  Spin  inversion.  For  an  open  shell  molecule,  the  code  should  give  the  same  results 
with  the  spins  of  all  electrons  flipped.  For  example,  a  triplet  state  should  give  the 
same  result  whether  the  two  unpaired  electrons  are  alpha  or  beta  spin. 

4.  Size  extensivity.  For  size  extensive  methods,  the  same  results  should  obtain  from 
a  single  molecule  and  from  two,  widely  separated  identical  molecules.  When  used 
with  very  tight  convergence  criteria  and  other  tolerances,  this  can  be  an  extremely 
accurate  test. 

5.  Separate  implementation.  After  all  other  tests  are  satisfied,  it  may  be  useful  to  make 
a  very  simple,  straightforward  implementation  of  the  equations  as  an  independent 
check.  For  example,  the  basic  equations  can  be  followed  literally  by  forming  a  large 
tensor  of  antisymmetrized  integrals  from  a  number  of  integral  lists,  and  likewise 
for  the  amplitudes,  and  comparing  the  results  element-by-element  with  the  primary 
implementation.  As  this  approach  tends  to  be  memory  intensive  and  slow,  small 
simple  test  cases  are  usually  the  best. 

There  are  also  a  number  of  tests  which  are  not  applicable  to  the  current  third-order 
triples  models,  but  should  be  kept  in  mind  for  situations  in  which  they  are  appropriate. 

1.  Active  space  invariance.  Fully  iterative  FSCC  methods  should  give  the  same  result 
for  any  choice  of  active  space. 
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2.    Invariance  to  orbital  rotations.  For  methods  which  (theoretically)  have  this  property, 
it  can  be  tested  by,  for  example,  using  an  ROHF  reference  function  with  and  without 
semicanonical  orbitals,  or  by  other  methods  which  modify  the  occupied-occupied  or 
virtual-virtual  parts  of  the  Fock  matrix. 
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